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This  dissertation  describes  a  model  for  acoustic  propagation  in  inhoihogeneous  fluid 
media,  and  explores  the  focusing  by  arrays  onto  targets  under  various  conditions.  The 
work  explores  the  use  of  arrays,  such  as  the  phase-conjugate  array,  for  underwater 
and  biomedical  applications.  Aspects  of  propagation  and  phasing  which  can  lead  to 
reduced  focusing  effectiveness  aje  described.  Among  the  most  important  debilitating 
effects  studied  here  are  medium  absorption,  medium  nonlinearity,  and  imperfect  initial 
phasing  of  the  signals  at  the  array  elements. 

Acoustic  wave  propagation  in  fluid  media  is  modeled  by  obtaining  a  wave  equation 
from  the  basic  equations  of  fluid  mechanics,  and  some  description  of  the  propagating 
environment  and  its  boundaries.  The  acoustic  wave  equation  couples  the  wave  motion 
to  the  medium’s  scalar  and  vector  time- varying  properties.  Bulk  velocity,  sound  speed, 
density,  attenuation  coefficient  and  nonlinearity  parameter  are  all  generally  functions 
of  three-dimensional  space  as  well  as  time  in  an  inhomogeneous  medium. 

The  present  study  uses  analysis  and  numerical  simulations  to  study  the  behavior 
of  the  acoustic  focusing  systems  described.  The  finite-difference  time-domain  (FDTD) 


method  is  used  to  solve  the  wave  equation  for  some  applications  in  underwater  and 
biomedical  acoustits.  The  nature  of  the  propagation  of  acoustic  disturbances  in  space 
and  time  in  a  continuous  medium  suggest  the  numerical  methods  used  to  solve  them 
here.  The  acoustic  pressure  disturbances  are  communicated  along  spatial  and  tem¬ 
poral  grids  in  a  natural  fashion.  The  strengths  and  weaknesses  of  the  FDTD  method 
are  discussed. 

Beyond  modeling  and  simulating  the  propagation  and  focusing  of  acoustic  fields, 
this  dissertation  looks  at  the  heating  effects  of  focused  ultrasound  in  an  absorbing 
thermoviscous  fluid.  The  application  considered  is  the  deposition  of  ultrasonic  en¬ 
ergy  onto  target  tissue  regions  with  the  purpose  of  affecting  therapeutic  heating  for 
hyperthermia.  The  acoustic  model  and  a  thermal  model  for  tissue  are  coupled  to 
solve  for  transient  and  steady  temperature  profiles  in  tissue-like  media.  Conclusions 
are  presented  on  the  effect  of  absorption,  nonlinearity,  and  temperature-dependent 
medium  properties. 
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Chapter  1 

INTRODUCTION 

The  focusing  of  intense  acoustic  beams  onto  targets  can  be  used  to  advantage  in 
military,  industrial,  and  biomedical  applications.  It  may  seem  ironic  that  the  utility 
of  the  intense  acoustic  pressures  generated  by  such  devices  is  often  in  its  useful  de¬ 
structive  abilities.  For  example,  focused  sound  may  provide  fast,  safe  minesweeping 
capability  for  ships  at  sea.  Focused  ultrasonic  fields  have  been  successfully  used  to 
break  off  oxidants  and  contaminants  from  semiconductor  surfaces  without  the  use  of 
solvents  or  direct  contact  by  abrasives.  Further,  the  heat  generated  near  regions  of 
focused  ultrasound  has  found  application  in  welding  plastic  parts  in  industry,  and  in 
treatment  of  deep-seated  tumors  in  tissue.  In  all  of  these  examples  focused  sound  was 
used  to  cause  physical  change  to  the  region  of  application,  most  commonly  causing 
permanent  destruction  of  the  material  at  the  focus.  This  aspect  is  shared  with  lasers 
in  that  intense  laser  beams  are  often  used  in  industrial,  medical,  and  military  appli¬ 
cations  to  destroy  objects  in  the  beam.  Contrast  this  usage  with  the  nondestructive 
diagnostic  and  imaging  uses  of  sound  and  laser  light  which  are  the  most  common 
usage  for  these  beams  at  low  intensities. 

This  study  explores  the  propagation  of  intense  sound  beams  from  focused  sources 
to  propagate  intense  beams  to  a  region  of  interest.  Applications  to  focusing  from  phase 
conjugate  arrays,  and  biomedical  use  for  therapeutic  tissue  heating  are  considered. 
The  emphasis  of  the  study  is  on  the  physics  of  the  wave  propagation,  the  debilitating 
effects  on  focusing  ability  which  may  arise,  and  the  modeling  and  behavior  of  tissue- 
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like  media  near  the  focus  of  ultrasonic  beams. 

In  Chapter  2  a  model  for  acoustic  wave  propagation  will  be  developed.  The  model 
is  used  to  study  the  motion  of  ultrasonic  pulses  propagating  in  inhomogeneous  fluid 
media.  First  we  consider  wave  motion  in  a  quiescent,  inhomogeneous,  lossy,  nonlin¬ 
ear  thermoviscous  fluid  medium.  Next,  propagation  in  non-quiescent  (time-varying) 
media  is  considered.  Since  the  propagation  can  be  very  complicated,  computer  sim¬ 
ulations  are  used  to  study  the  behavior  of  the  acoustic  field,  especially  near  the  focal 
zone.  The  finite-difference  time-domain  (FDTD)  method  for  solving  the  wave  equa¬ 
tion  is  briefly  outlined  at  the  end  of  the  chapter  to  facilitate  the  presentation  of  some 
numerical  examples  from  the  model. 

Chapter  3  presents  formal  derivations  for  the  field  due  to  point  sources  using  the 
Green’s  function.  The  analysis  is  extended  to  arrays  of  small  sources  and  extended 
sources  in  free  space.  The  superposition  of  linear  acoustic  fields  is  used  to  illustrate 
the  need  for  geometric  focusing  and  phasing  to  form  intense  foci.  A  discussion  of 
beamsteering  and  some  common  configurations  for  focused  acoustic  sources  are  pre¬ 
sented  as  well. 

Chapter  4  considers  time  reversal  or  phase  conjugation  arrays.  Such  arrays  have 
been  shown  to  have  remarkable  abilities  to  automatically  compensate  for  the  effects 
of  complicated  propagation  paths,  scattering,  and  multipath  effects  in  fluid  channels. 
In  this  study  we  demonstrate  analytically  and  through  numerical  experiments  that 
the  robustness  of  time  reversal  arrays  for  retrofocusing  holds  even  under  nonlinear 
propagation  conditions.  The  caveat  to  this  last  statement  is  that  the  propagation 
medium  must  be  nearly  lossless,  and  no  modification  to  the  signal  other  than  time 
reversal  may  be  applied  at  the  array.  The  effect  of  various  debilitating  factors  on  the 
ability  of  such  arrays  to  form  an  intense  focus  at  a  target  location  is  studied.  For 
high-intensity  purposes,  absorption,  amplification  at  the  array  elements,  and  error 
in  phasing  will  compromise  the  time  reversal  invariance  necessary  for  effective  time 
reversal  focusing.  All  of  these  effects  act  to  degrade  the  focusing  of  time  reversal 
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systems,  and  are  shown  to  be  especially  detrimental  for  finite-amplitude  propagation. 
An  effort  is  made  to  elucidate  the  effects  of  individual  propagation  conditions  on  the 
propagation  and  array  focusing,  so  minimizing  the  number  of  simultaneously  acting 
phenomena  will  minimize  confusion  as  to  cause  and  effect.  So  for  example,  in  studies 
investigating  the  effect  of  phase  error  on  the  focusing  performance,  the  absorption 
and  nonlinearity  are  suppressed  so  as  not  to  cloud  the  issues  being  studied  in  that 
example. 

In  Chapter  5  the  effects  of  focused  ultrasound  in  tissue  are  studied.  A  thermal 
model  for  tissue-like  media,  based  on  the  bioheat  equation,  is  coupled  to  the  pre¬ 
viously  described  acoustic  model.  A  heat  equation  is  solved  where  the  deposition 
of  thermal  energy  into  the  tissue  from  focused  ultrasound  is  countered  by  conduc¬ 
tion  and  perfusion  heat  losses.  It  is  known  that  material  properties  such  as  sound 
speed  and  absorption  coefficient  in  biological  media  are  temperature-dependent.  The 
effect  of  the  acoustic  heating  on  the  thermal  and  acoustic  properties  of  tissue  will 
have  a  feedback  effect  on  the  sound  field  and  rate  of  heat  deposition.  The  effect  of 
time/temperature- varying  sound  speed  and  absorption  coefficient  is  studied  in  sim¬ 
ulations,  using  published  laboratory  data  for  these  parameters  as  a  function  of  tem¬ 
perature.  It  is  found  that  the  sharp  increase  in  absorption  coefficient  as  a  heating 
treatment  progresses  will  have  an  important  effect  on  the  rest  of  the  treatment  by 
accelerating  the  heating  near  the  focal  zone. 

Each  chapter  is  individually  summarized  at  its  end,  but  Chapter  6  ties  the  disser¬ 
tation  together  by  presenting  conclusions  from  the  findings  of  the  whole  study  in  a 
coherent  fashion.  Chapter  6  also  describes  some  outstanding  questions  not  addressed 
in  this  dissertation,  and  proposes  some  directions  for  extending  this  work.  Some  of 
the  most  important  points  not  covered  in  the  present  study  are  the  role  of  moving 
medium  in  time  reversal  systems,  the  role  of  time- varying  properties  of  tissue  other 
than  sound  speed  and  attenuation  (for  example,  density,  nonlinearity  coefficient,  vis¬ 
cosity),  and  the  onset  of  phase  change  to  treated  tissue  following  the  denaturing,  or 
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cavitation  inception.  Only  the  phase  errors  in  the  initial  phasing  of  time  reversal 
systems  were  studied.  Phase  errors  in  the  waveforms  themselves  can  also  be  expected 
in  real  systems.  Given  the  present  interest  in  the  ability  of  forming  high-intensity  foci 
using  time  reversal  arrays,  the  time  reversal  systems  are  evaluated  here  based  on  the 
intensity  of  the  focal  spot  generated,  not  on  the  fidelity  of  reproduction  of  the  initial 
waveforms.  This  may  not  be  the  preferred  judgement  criterion  in  some  of  the  “softer” 
uses  of  time  reversal  arrays,  such  as  communication  or  imaging. 

The  present  study  relies  heavily  on  numerical  simulations.  While  all  the  algo¬ 
rithms  and  computer  codes  used  in  this  study  are  original,  the  numerical  methods 
used  are  not  new.  Hence  the  description  of  the  numerical  techniques  used  is  dis¬ 
cussed  in  the  Appendix  rather  than  in  the  body  of  the  dissertation.  Some  comments 
regarding  the  simulation  graphical  output  are  due:  to  simplify  reproduction  of  this 
dissertation,  gray-scale  coloring  was  used  for  all  figures.  Unfortunately,  this  compro¬ 
mised  some  of  the  detail  visible  in  the  original  color  output.  A  few  colored  versions 
of  this  dissertation  were  produced,  but  the  official  document  remains  in  black  and 
white.  Also,  some  insight  was  gained  by  generating  and  viewing  animated  sequences 
of  acoustic  pulse  propagation,  especially  for  time  reversal  and  therapeutic  ultrasound 
pulse  propagation.  These  too  remain  axchived  with  the  author  and  could  not  be  in¬ 
cluded  in  the  official  dissertation.  The  simulations  presented  were  carried  out  on  the 
Silicon  Graphics  (Origin  2000)  parallel  computer  at  Boston  University  using  the  Cray 
parallellizing  FORTRAN-77  compiler.  Postprocessing  was  done  on  the  Origin  2000 
and  Apple  Macintosh  computers.  No  Intel  Pentium  chips  were  used  in  this  study. 
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Chapter  2 

MODELLING  THE  ACOUSTIC  WAVE 

The  first  step  to  simulating  an  acoustic  field  is  to  decide  on  a  model  for  sound 
propagation.  This  chapter  derives  the  model  wave  equation  from  the  basic  equations 
of  fluid  mechanics  and  thermodynamics  for  a  thermoviscous  fluid.  This  study  concerns 
itself  with  the  compressional  acoustic  wave  action  in  fluid  media.  The  existence  or 
generation  of  shear  waves  is  not  considered.  The  wave  equation  of  acoustics  can  be 
derived  from  the  fundamental  equations  of  fluid  mechanics  and  thermodynamics.  In 
order  to  retain  the  nonlinear  behaviour  of  the  wave,  it  is  necessary  to  retain  the 
nonlinear  terms  which  arise  in  the  derivations.  We  keep  terms  up  to  second  order  in 
this  study,  a  regime  referred  to  as  finite-amplitude  acoustics.  Keeping  terms  of  higher 
than  second  order  results  in  much  more  complicated  expressions. 

2.1  The  Basic  Equations  of  Fluid  Mechanics  and  State  in  Quiescent 
Fluids 

The  equations  used  to  derive  the  wave  equation  in  fluids  are  Euler’s  equations  (con¬ 
tinuity  and  momentum),  and  an  equation  of  state  relating  pressure  to  density  and 
entropy.  The  derivations  in  this  section  apply  to  quiescent  fluids^,  the  case  of  fluids 
with  time-varying  background  properties  will  be  derived  later  in  this  chapter.  The 
continuity  equation  expresses  the  conservation  of  mass  in  a  fluid  volume,  and  is  given 


^Defined  by  Pierce  [61]  p.l4  as  a  fluid  whose  background  properties  are  not  dependent  on  time, 
and  whose  background  velocity,  uo,  is  zero. 
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by  Kundu  [51]  as 

|^  +  V-(pu)  =  0,  (2.1) 

where  the  density,  p,  and  the  velocity,  u,  are  both  functions  of  space  and  time. 
Acoustic  disturbances  are  considered  to  be  small  perturbations  to  the  background 
state  of  the  fluid,  and  the  density,  velocity,  and  pressure  are  written 


p 

=  Po  +  p'  1 

(2.2) 

u 

=  Uo  + 

(2.3) 

p 

=  Po  +  p'- 

(2.4) 

The  quantities  with  the  subscript  zero  refer  to  the  fluid’s  background  (ambient)  prop¬ 
erties  and  are  only  functions  of  space,  while  the  quantities  having  the  superscript 
prime  refer  to  perturbations  from  the  background  state  due  to  the  acoustic  wave,  and 
are  generally  functions  of  both  space  and  time.  Since  Uq  =  0  in  all  of  the  cases  treated 
in  this  study,  u  =  u'. 

Expanding  (2.1)  and  using  (2.2)  we  obtain  the  continuity  equation  for  a  quiescent 
inhomogeneous  fluid, 

^  +  poV  •  u'  =  -p'V  •  u'  -  u'  •  Vp'  -  u'  •  Vpo.  (2.5) 

iJv 

In  a  homogeneous  medium,  the  last  term  in  (2.5)  would  be  dropped. 

The  next  piece  of  fluid  mechanics  we  will  use  to  derive  the  wave  equation  is  the 
momentum  equation, 

P^  +  Vp-Fb-{\  +  2p)  V(V  •  u)  =  0,  (2.6) 

where  p'  is  the  acoustic  pressure,  t  is  time,  and  p  and  A  are  the  viscous  Lame  coeffi¬ 
cients.  D  denotes  a  material  derivative. 


D 

Dt 


V. 


(2.7) 
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The  body  force  Fb  and  the  gradient  of  the  background  pressure,  Vpo  cancel,  since  in 
the  absence  of  sound,  p'  =  u  =  0,  and  so  (2.6)  is  merely 

Vpo  =  Fb>  (2-8) 

A  question  may  arise  as  to  the  importance  of  each  term  in  the  equations  of  fluid 
mechanics  to  the  overall  behavior  of  the  fluid.  To  answer  this  question  and  to  elucidate 
more  of  the  nature  of  the  acoustics  we  use  an  ordering  scheme  that  clarifies  which 
terms  are  small  compared  to  the  others.  In  this  study  terms  of  the  order  three  and 
higher  are  ignored. 

2.2  On  the  Ordering  of  Terms 

We  shall  describe  the  basis  for  ordering  the  terms  encountered  throughout  our  analy¬ 
sis.  We  assume  that  background  parameters  (e.g.  po,  po  Co)  are  of  order  one,  denoted 
by  0(1).  Disturbances  due  to  the  acoustic  waves  are  considered  to  be  order  epsilon, 
0(e),  and  are  much  smaller  than  the  background  values  and  0(1)  terms.  These 
include  pf,  u',  p' ,  etc.  The  viscous  (Lame)  coefficients  are  also  considered  small  com¬ 
pared  to  the  0(1)  terms,  and  are  first  order  terms  denoted  by  0(t]).  Differential 
changes  in  background  parameters  (e.g.  ^  and  Vpo)  are  assumed  to  be  C>(C)  per¬ 
turbations,  where  e,  rj  and  (  are  all  similarly-small  and  much  less  than  one.  Table  1.1 
summarizes  the  ordering  scheme  used  in  this  study.  Note  that  the  product  of  two  first 
order  terms  is  a  second  order  term,  thus  allowing  the  substitution  of  first  order  terms 
into  relations  that  are  second  order  while  maintaining  second  order  accuracy.  Note 
the  distinction  between  the  terms  “order  one”,  0(1),  and  “first  order”,  0(e).  Order 
one  terms  could  be  called  “zeroth  order”  terms.  Some  examples  using  the  ordering 
scheme  described  are 
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Table  2.1:  The  ordering  scheme  used  in  this  study  (n  >  0,  m  >  1).  e,ri  and  (  are 
much  smaller  than  1.  s  represents  space  or  time  in  the  differentials. 


^n^/m 

ds^ 

8s" 

8s" 

8s" 

8s" 

0(e™) 

C>(e“) 

0{c^) 

o(C”) 

0(C") 

^(^n+m) 

0(77"+™) 

— =  (9(1)  0(1“^)  =  (9(1),  (zeroth  order) 

PoCo 

=  0(1®)  0(e)  =  0(e),  (first  order) 

—  0(1^)  0{€^)  =  C*(e^),  (second  order) 

=  C>(e^)  (9(C)  =  0{e^  C)  (third  order). 

Manipulation  of  second  order  relations  such  as  (2.5)  is  possible  using  the  first  order 
relations  described  in  [13,  68]  while  maintaining  second  order  accuracy.  For  our  pur¬ 
poses  we  use  the  first  order  linearized  relations  for  mass  and  momentum  conservation. 


^  (2.10) 
=  -Vy,  (2.11) 


as  well  as  a  first  order  relation  describing  the  dependence  of  density  fluctuation  on 
the  background  sound  speed  and  the  acoustic  pressure. 


,  1  / 

p  =:aP- 


(2.12) 


These  relations  allow  us  to  rewrite  the  right  hand  side  of  (2.5)  to  obtain  a  more  useful 
form  of  the  continuity  equation. 


dp'  ^  ,  1  dp^ 

dt  pqCq  dt 

which  will  be  used  to  derive  the  acoustic  wave  equation. 


(2.13) 


9 


2.3  The  Absorbing  Nonlinear  Wave  Equation 


A  detailed  treatment  of  the  wave  equation  for  nonlinear  acoustics  in  fluids  can  be 
found  in  Hamilton  and  Blackstock  [42].  To  obtain  the  acoustic  wave  equation,  the 
time  derivative  of  equation  (2.13), 


./2 


dy  „  1  ay 

dt^  at  poci  at^ 


is  subtracted  from  the  divergence  of  equation  (2.6) 


_  an' 

■  aT’ 


(2.14) 


PoV  •  ^  ^  •  V/r)o  +  VV  =  (A  +  2ix)V\V  •  u')  +  V(V  •  u') .  V(A  +  2p), 


to  yield 


(2.15) 


Qw! 

V/)o--^  =  0.  (2.16) 


We  return  to  the  first  order  relations  (2.10)  and  (2.11),  and  instead  of  (2.12)  we  use 
the  state  equation  and  entropy  equations  to  second  order,  derived  in  Hamilton  and 
Blackstock  [42] 


.  P 

P=^ 


K 


1 


'  W  -  ~Ap’\ 


^0  po^o  ^v)  PoCqSA 


(2.17) 


where  k  is  the  thermal  conductivity.  The  constants  A  and  B  are  the  coefficients  of 
the  polynomial  describing  the  relationship  between  the  fluid’s  density  and  acoustic 
pressure  variation  described  by  Beyer  [6].  B  j A  is  proportional  to  the  ratio  of  coeffi¬ 
cients  of  the  quadratic  and  linear  terms  in  the  Taylor  series  expression  for  pressure 
as  a  function  of  condensation.  We  use  the  second  time  derivative  of  (2.17), 

a^p'  1  ay 


at^ 


4  at^ 


K 

Pocl 


/I  i\a^p'  B  ay^ 
at^  2Apo4  at^ ' 


(2.18) 


to  eliminate  acoustic  density  from  the  nascent  wave  equation  (2.16).  We  are  now 
left  with  an  equation  with  the  acoustic  pressure,  p',  as  the  sole  dependent  variable. 
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Collecting  like  terms,  (2.16)  can  be  written  as 


- 2  ^l2 - ^ - 2  (  7^ - 7^  )  “Sis’ 

Cq  dt  po  Po^O  \(^v  ^pj 


B  dY  ,  1 

2ApoCq  dB  ^  ^  PoCq  dt^  poc^  dt"^ 


(2.19) 


The  terms  in  the  above  equation  are:  (a)  the  D’Alembertian,  present  in  all  wave 
equations.  This  term  describes  the  propagation  of  a  pulse  in  time  and  space;  (b) 
the  ambient  inhomogeneity  in  the  medium’s  density,  and  is  zero  for  a  homogeneous 
medium;  (c)  and  (e)  the  loss  terms,  due  to  the  thermal  conduction  and  the  viscosity 
of  the  fluid  respectively;  and  finally  (d)  and  (f),  the  nonlinear  terms  arising  from  the 
equation  of  state  and  the  continuity  equation  respectively.  The  nonlinear  terms  can 
be  combined  using  the  nonlinearity  coefficient,  /?,  which  is  related  to  the  number  Bf  A 


^  =  1  + 


(2.20) 


Simplification  allows  us  to  write  the  wave  equation  in  terms  of  the  acoustic  diffusiv- 
ity,  5,  which  accounts  for  both  thermal  and  viscous  losses,  in  a  form  attributed  to 
Westervelt  [73], 

W  -  +  -4^  =  0-  (2-21) 

eg  av  po  cl  ot^  pqcI  dv 

We  now  have  a  second-order  wave  equation  describing  the  acoustic  pressure  in  terms 
of  space,  time,  and  the  fluid’s  material  properties. 

The  simulations  presented  in  this  study  use  an  absorption  coefficient,  a,  which  is 
related  to  the  acoustic  dilfusivity  by 


and  has  the  units  of  Nepers/meter  [Np/mj. 


(2.22) 
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2.4  Similarity  Solution  for  the  Wave  Equation  in  Quiescent  Fluids 


The  wave  equation  (2.21)  is  given  in  its  full  dimensional  form.  Actual  values  for 
Co,  ao)  Poi  0Oi  and  source  pressure,  Psource;  must  be  provided  for  the  specific  prop¬ 
agation  medium  under  consideration  to  obtain  a  solution  for  that  case.  A  useful 
technique  that  can  be  used  to  gain  insight  into  the  essential  nature  of  the  wave  equa¬ 
tion  is  a  similarity  analysis.  The  advantage  of  using  similarity  analysis  is  that  two 
systems  with  different  background  parameters  and  driving  pressures  can  be  studied 
using  the  same  solution,  provided  that  certain  nondimensional  quantities  that  char¬ 
acterize  the  wave  equation  are  the  same  for  the  two  systems.  This  reduction  of  a 
differential  equation  to  its  simplest  form  with  no  explicit  dimensional  variables  is  also 
known  as  nondimensionalizing  the  equation.  In  this  section  we  nondimensionalize  a 
simplified  version  of  (2.21)  for  illustration  purposes  in  polar  cylindrical  coordinates 
in  a  homogeneous  medium.  The  dimensional  form  of  the  wave  equation  is 


V  V  -  -2 


1  2a  av 


(2.23) 


A  characteristic  set  of  nondimensional  variables  with  hats  [e.g.  i)  is  used  instead 
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of  the  dimensional  variables  {e.g.  t)  used  in  (2.23). 

.  P' 

7, - ’ 

/^source 
A  ‘f' 

Ao 

a:  =  p 

i  =  tu, 


c-  — 
Co 


(2.24) 


P  = 


a  = 


Pq 

a 

cto' 


The  zero  subscripts,  {e.g.  po)  denote  background  values.  In  this  case  we  take  the 
background  parameters  to  be  constants,  as  the  medium  is  assumed  to  be  spatially 
homogeneous  and  quiescent. 

In  polar  cylindrical  coordinates  the  scalar  Laplacian  is 


1  d  .  dp'  5  V  _  1  ^P'  ^  V 

r  dr  dr  ^  dx'^  dr"^  r  dr  dx'^  ’ 


(2.25) 


where  we  have  assumed  axial  symmetry  exists  and  neglected  any  derivatives  in  az¬ 
imuthal  angle,  6.  This  is  a  useful  geometry  found  in  most  experimental  bowl  trans¬ 
ducers  in  research  laboratories.  The  operators  in  (2.23)  are  nondimensionalized  using 
the  angular  wavelength 


Ao 


Co 

? 

(Jj 


(2.26) 
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and  the  transformations  on  the  operators 


dr  d.r  1 

dr  dr  Aq  Aq’ 

dp  _  Psource  dp 

dr  Ao  dr' 

dy  _  p'  d'^p 

dr'^  Aq  dr'^ ' 

dp  _  Psource  dp 

dx  Aq  dx ' 

dy  _  p'  d'^p 

dx'^  Ag  ’ 

1  1 


di 

Ft 


Aor’ 

d,  . 

=  jfiu,)  =  a,, 


dp' 


dp 


dt 


dy 

dt^ 

d^p' 

dt^ 


2  ^ 
^  ^source  ? 

Ot^ 


—  Ct7  ^source 


di^' 


(2.27) 


Applying  the  above  transformations  and  multiplying  through  by  gives 

a  dimensionless  form  of  the  wave  equation  in  a  homogeneous  medium  in  cylindrical 
coordinates 


d'^p 

dP 


1  dp  d'^p  ^^P  ,  o  \  I  2  ^0  Psource 

f  dr  dx"^  dP  °  °  dP  Po  Cq 


We  pause  now  to  examine  the  remaining  factors  in  front  of  our  operators  in  the 
nondimensional  wave  equation  (2.28).  Two  points  of  primary  importance  are  noted: 


1 .  Two  nondimensional  terms  are  all  that  controls  the  behavior  of  any  particular 
case  of  this  equation,  so  long  as  the  nondimensional  factors  are  the  same.  These 
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factors  are 


ao  A 


0) 


and 


(2.29) 


Po  Psource 

PQ  Co 


—  Po  C- 


Note  that  the  absorption  term  contains  a  factor  called  the  absorption  length 


otQ  Ao,  which  is  the  amount  of  absorption  suffered  in  propagating  the  distance 
Ao-  Also  note  that  the  acoustic  Mach  number, 


^source  _  ^source 

6  —  ^ , 

Co  PO  Cq 

appears  as  a  multiplier  in  front  of  the  nonlinear  term. 


(2.30) 


2.  The  ratio  of  the  two  nondimensional  parameters  in  (2.29)  is  a  measure  of  the 
importance  of  nonlinearity  to  absorption.  This  ratio  is  derived  by  Blackstock 
[7],  and  is  called  the  Gol’berg  number,  F, 

r  =  (2.31) 

oioAo 

In  the  coming  chapters  we  will  see  how  the  nonlinearity  and  the  absorption  both 
come  into  play  in  the  propagation  of  finite-amplitude  waves,  and  how  the  nonlinear 
steepening  and  the  absorption  are  competing  mechanisms  in  nonlinear  wave  propaga¬ 
tion.  Blackstock  points  out  that  the  Gol’berg  number  can  be  used  as  an  estimate  of 
the  overall  nature  of  a  pulse’s  propagation.  If  F  >  1  then  nonlinearity  dominates  and 
the  waveforms  will  tend  to  develop  steep  shocked  profiles,  if  F  1  then  absorption 
dominates  and  the  waveform  decays  without  severe  steepening. 


2.5  The  Wave  Equation  with  Time-Varying  Background  (TVB)  Param¬ 
eters 

We  re-derive  the  wave  equation  from  the  nonlinear  equations  of  fluid  mechanics  and 
state,  now  allowing  the  background  sound  speed  and  density,  po{t)  and  co{t),  to  be 
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time-varying.  We  refer  to  this  as  the  time-varying  background  (TVB)  scenario  or 
solution.  The  motivation  for  deriving  the  TVB  wave  equation  is  that  changes  which 
occur  to  the  background  values  of  po  and  cq  at  time  scales  comparable  to  those  of 
the  acoustic  variables  may  contribute  to  the  wave  equation.  The  fluid’s  density  and 
nonlinearity  coefficients  axe  assumed  to  remain  constant.  This  constraint  is  born  of 
the  fact  that  only  the  primary  propagation  variable  (sound  speed)  and  the  primary 
heating  parameter  (absorption  coefficient)  have  been  reported  for  tissue  as  a  function 
of  temperature.  Thus  for  this  analysis  the  density  and  nonlineaxity  coefficient  are 
assumed  to  be  constants. 

The  continuity  equation  (2.13)  and  the  momentum  equation  (2.6)  are  manipulated 
as  in  Section  2.3,  except  that  time-derivatives  of  po{t)  and  Co(t)  axe  assumed  to  be 
non-zero.  These  derivatives  are  taken  to  be  first  order  terms.  The  continuity  equation 
(2.1)  is  now 

%  +  ^  +  PoV  •  u  =  -p'V  •  u  -  u  •  Vp'  -  u  •  Vpo.  (2.32) 

at  at 

Equation  (2.32)  to  first-order  accuracy  (neglecting  0(e^)  and  higher-order  terms)  is 

^  +  ^  +  /-oV.u  =  0  (2.33) 


from  which  we  derive  the  first-order  relations 


V-u  = 

Po  dt  Po  dt 


Next  we  consider  the  momentum  equation, 


p^  +  Vp  -  Ffl  -  (A  -I-  2p) V(V  •  u)  =  0. 


(2.34) 

(2.35) 

(2.36) 


In  the  absence  of  sound,  p'  =  u  =  0,  and  so  the  gradient  of  the  background  pressure 
and  the  body  force  cancel  as  in  (2.8).  So  (2.36)  for  the  TVB  case  to  second  order 
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accuracy  is 


Po-^  +  W  =  -p'^  -  pou  •  Vu  +  (A  +  2ju)v^u, 


dt  ‘  ^  dt 

where  we  have  used  the  vector  identity 


V(V  •  u)  =  V^u  +  V  X  V  X  u 


(2.37) 


(2.38) 


and  assumed  that  the  flow  is  irrotational,  vis.  V  x  u  =  0.  The  momemtum  equation 
to  first-order  accuracy  in  a  TVB  fluid  can  be  written  to  first  order 


at  po 


(2.39) 


Note  that  this  is  identical  to  that  of  the  quiescent  medium. 

A  first  order  wave  equation  can  be  obtained  by  subtracting  the  time-derivative  of 
(2.33)  from  the  divergence  of  (2.39), 


vy- 


dt^ 


=  0. 


The  equation  of  state  for  a  TVB  fluid  is 
P‘ 

To  first  order  this  is 


(2.40) 


(2.41) 


(2.42) 


p'  =  cIp'  +  0(6^), 

from  which  the  first-order  relations  below  may  be  derived 

p'  =  ip'  +  0{^), 

dt  eg  af  ' 

Vp' =  ivp' +  C'Ce''). 

% 

The  next  step  towards  obtaining  a  second-order  wave  equation  for  nonlinear  ab¬ 
sorbing  TVB  media  is  to  use  first-order  relations  to  eliminate  p'  and  u  from  the  wave 


(2.43) 
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equation,  and  to  express  the  equation  in  terms  of  p' .  We  write  the  continuity  equation 
(2.32)  as 


dpo  ,  dp\  ^  p'  dpo  ^  p'  dp'  , 

^  +  ^0^  ■  “  =  -  u  •  Vp  -  u  •  Vpo  -  u  •  Vpo 

at  at  po  at  po  at 


p'  dpo 


+ 


pocl  dt  2poCo  dt 


1  dp'^  u  , 

-  3  •  Vp  -  u  •  Vpo 


P'  9po  ^  1 


dp‘ 


/2 


+ 


c5 

Po  du^ 
2co  dt 


(2.44) 


u  •  Vpo  "H  ^(^^)- 


PoCq  dt  2poCo 

If  we  define  the  Lagrangian  density  (c./.  Hamilton  and  Blackstock  [42]  Chapter  3) 


1  2 

2'^“  2p„c§’ 


(2.46) 


then  the  quantity 


Cq  dt 


Po  dv?  1  5p' 


./2 


2co 


2poCo  dt 

so  the  continuity  equation  can  be  written  as 


+  0{e% 


(2.46) 


dpo  ,  dp'  ^  „  _  p'  dpo  ,  ^  ,  1 

dt  dt  PoCq  dt  dt  Cq  dt 


(2.47) 


Equation  (2.47)  is  the  same  as  that  would  be  obtained  in  a  quiescent  fluid  with  the 
exception  of  the  terms  containing  which  are  due  to  the  TVB  character  of  po- 
We  also  rewrite  the  momentum  equation  for  a  TVB  fluid  (2.36)  using  the  0(e) 
substitutions  (2.35),  (2.39),  and  (2.43), 


dvL  „  ,  ^  1  dpo 

_  A  +  2p  ap' 
dt'^ 


i^)+ 

Po  dt  pocl 


1 


Po 


2/^0  ^0 

A  +  2p  ap' 

Po  at 


Vp' 


/2  Po 


(2.48) 


where  the  gradient  of  the  Lagrangian  density  (2.45)  is 

V£  =  ^ 

2  2^(,cJ 


,2  • 


(2.49) 


18 


To  obtain  a  wave  equation  we  must  subtract  the  time-derivative  of  (2.47), 


d'^p'  d 


from  the  divergence  of  (2.48), 


^Vo 


Vu 


dpo 


df^  '  ”  dt 

1  dpo  dp'  1  aV*  1  d^c 

A  ,o  4” 


(2.50) 


PoCq  dt  dt  pqc^  dt^ 


4  dt'^ 


vv+p„|v.u=-(4±Mv^|:-v^£, 

dt  po4  ot 


(2.51) 


to  obtain,  after  substitutions  and  algebra, 
d^p' 


vV 


dt^ 


+ 


1 

d^p'^ 

(A  +  2/i) 

po4 

dB 

1 

o 

d'^po 

,_1| 

(dpoV 

dt^ 

Po 

[dtj 

dt  po 
2  dp^ 


po4 


I  dt  V 


dt'^ 


(2.52) 


TVB, 


PO 


This  is  the  “nascent”  wave  equation  for  the  TVB  case.  The  terms  in  (2.52)  are  similar 
to  those  derived  for  the  quiescent  case,  with  the  addition  of  the  terms  labeled  TVBp^ 
which  are  due  to  the  TVB  nature  of  po-  The  background  sound  speed’s  contribution 
to  the  TVB  terms  will  come  about  in  the  following  equations  where  the  equation  of 
state  is  used  to  eliminate  the  p'  from  (2.52). 

Applying  first-order  substitutions  (2.43)  to  the  state  equation  (2.41)  gives  a  second- 
order  expression  for  p'  as 
.  1 


iP 


n  ,  ..si  dp  1  B  ,2  . 


c, 


Po% 


■^2A 


p'^  +  0{e^). 


(2.53) 


Differentiation  with  respect  to  time  twice  gives 

1  ay  1  B  d^p'^  1  «(7-1)  5V 

~dF~4dt'^  Po42A  dB  po4  Cp  dt^  4  dt  dt  ' 

TVBco 


(2.54) 


where  the  term  labeled  TV Bcq  is  a  result  of  the  TVB  nature  of  cq. 
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Substituting  (2.54)  into  (2.52)  elinainates  the  p'  variable  from  the  wave  equation, 
leaving  the  wave  equation  for  the  TVB  medium  as 


^  _  (A  +  2)»)  ay _ 1_K 

cl  dt'^  PoCq  V  2aJ  po(^  dt^  poc^ 


(7  -  1) 

Cp  dt^ 


+  TVB,,  +  TVB^ 


or  after  some  manipulation,  and  neglecting  the  C  terms. 


c4  Q^2  p^ 


Vp'  ■  Vpo 


d'^Po  ,}_( AA\  a.  j.  _  A 

dB  po\dt)  poCq  dt  dt  dt  dt 


(2.55) 


(2.56) 


Note  that  S  =  b/po  where  6  =  A  +  2//  +  ~  !)• 

As  a  comment,  truncating  higher  than  second  order  terms  in  the  derivation  of  the 
above  equations  results  in  a  reajsonably  compact  wave  equation.  Keeping  higher  order 
terms  than  second  order  would  result  in  significantly  more  complicated  equations. 

2.5.1  Nondimensionalizing  the  TVB  Wave  Equation 

We  nondimensionalize  (2.56)  for  a  one-dimensional  case  using  characteristic  nondi- 
mensional  variables  with  hats  {e.g.  po)  as  before,  with  some  new  nondimensional 
variables 


Po  =  — , 
Poo 

ao 

oo  = - , 

aoo 

/?oo' 


(2.67) 
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The  parameters  with  zero  subscripts,  (e.^.  po)  denote  background  values,  which  are 
functions  of  time^.  The  parameters  with  double-zero  subscripts  {e.g.  poo)  are  back¬ 
ground  parameters  of  the  undisturbed  background  state  of  the  fluid  before  any  inter¬ 
actions  altering  the  background  properties  take  place  in  the  fluid,  and  are  constants. 

The  wave  equation  in  nondimensional  form  is 

d'^p  1  d'^p  ^  do  \  ,  ^0  /5ooPsource  Op  dpo 

_  +  —f  ^  Q 

dP  Po\diJ  PoCq  di  di  Cq  di  di 
2.5.2  Ordering  the  Terms  in  the  TVB  Wave  Equation 

In  subsequent  chapters  on  thermal  effects  of  focused  ultrasound  we  will  find  it  conve¬ 
nient  to  show  that  not  all  terms  from  (2.58)  need  to  be  computed  if  some  are  much 
smaller  than  the  others.  To  this  end,  we  now  examine  (2.58) ’s  terms  to  set  their 
magnitudes  to  some  scale  for  comparison. 

For  reasonable  amounts  of  variation,  the  background  parameters  are  0(1), 

Co,  po,  do,  /3o  =  ^(1)5  (2.59) 


because  for  any  background  parameter  xo 


Xo  _  g(l) 

Xoo  0(1)' 


(2.60) 


We  can  thus  drop  the  hatted  groups  of  parameters  leading  the  terms  in  (2.58)  as  0{1), 


so 


S^p  d^p  ~  d^p  ,  /dooPsource  d'^p^  dp  dpo 

_  / ^PoV  dpdpo  ^ 

dP  ^  \di )  di  di  di  di 


^We  assume  that  some  mechanism  changes  the  background  properties.  In  the  case  of  focused 
ultrasound  surgery,  background  parameters  have  been  measured  to  be  functions  of  temperature. 
Since  the  temperature  of  the  tissue  in  this  case  is  time-dependent,  so  are  the  background  properties. 
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In  order  to  determine  the  importance  of  the  terms  relative  to  one  another,  we 
rewrite  (2.61)  as 

d'^p  d'^p  ,  X  ,  f^ooPsouvced'^p^  dpdpQ 


&^p  d'^p  \  ^^P  ,  /^OOPsource  d'^p^ 

dil  ^°dil  poocgo  dil 


dx  dx 


Tg  d'^po  (  Tgdpo 

ri  dil  V6  dib 


n  dig  dib  n  dig  dib  ~  ’ 


(2.62) 


where  two  nondimensional  times  have  been  defined  corresponding  to  the  acoustic  and 
the  thermal  time  scales: 


tg  —  , 

Ta 


(2.63) 

(2.64) 


and  we  used  the  manipulation  (an  example  is  given  here), 

dcp  _  d^^  _  dcp  d{tlTb)  _  d^T^ 

dig  dib  dig  dibd{tlTg)  dUn 

Now  that  the  equation  is  written  in  terms  of  partial  derivatives  which  are  all  order 
unity,  we  can  now  use  the  factor  Tg/rb  as  an  estimator  of  relative  magnitudes  of  the 


terms. 


0(1)  0(1) 


O(aooAoo) 


O  (/JooPsource  /  POO  Cqo  ) 


dp  d^p  ~d^p  /?00Psourc 

- +  ttoo^oo  TTfT  d - 2~ 

dx^  dtl  dtl  poocgo 


0(to/t6)^ 


0{TalTbf 


<^(Ta/T6) 


(2.66) 


0{TalTb) 


T^d^po  ,  (Tgdp 


dil 


n  dib 


Tb  dig  dib  n  dig  dib 


As  an  example,  we  perform  an  order-of-magnitude  analysis  of  (2.66)  for  the  case  of 
soft  tissue  medium  insonated  by  a  sinusoid  at  acoustic  angular  frequency  Ug  of  1  MHz. 
Our  experience  [37],  and  laboratory  measurements  [45,  11,  65]  at  typical  therapeutic 
biomedical  device  source  pressures  (1  MPa)  shows  that  for  an  acoustic  source  with 
period  Tg  =  Ips,  tissue  heating  occurs  over  time  scales  on  the  order  of  rj  =  Is.  Thus 


-  «  10-®. 

n 


(2.67) 
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For  typical  background  reference  properties  of  soft  tissue, 


Coo  =  1600  m/s, 

=  1100  kg/m^ 
ooo  =  4.5  Np/m, 
l3oo  =  5.5, 

we  conclude  that  the  D’Alembertian  terms  in  (2.66)  are  of  magnitude  unity,  the 
absorption  and  the  nonlinear  terms  are  of  magnitude  10“^,  and  the  TVB  terms  are 
of  magnitude  lO'®  and  10“^^.  This  result  shows  formally  that  for  hyperthermia 
applications  (to  be  discussed  in  detail  in  Chapter  5)  the  TVB  is  negligible  in  its  effect 
on  the  wave  equation.  In  other  words,  the  acoustic  field  is  not  measurably  influenced 
by  the  change  in  tissue  background  parameters  when  viewed  at  time  scales  Ta. 


2.5.3  Conclusions  Regarding  the  Wave  Equation  for  TVB  media 

The  analysis  suggests  that  dynamic  background  medium  properties,  po{t)  and  co(it), 
can  be  important  enough  to  be  included  in  a  second-order  wave  equation.  The  nondi- 
mensionalized  equations  suggest  that  the  time-variation  of  po  and  cq  can  be  neglected 
if  the  characteristic  time  of  change  for  these  parameters  is  much  longer  than  the 
acoustic  characteristic  time  scale,  the  period.  From  this  we  conclude  that  while  in 
the  biomedical  ultrasound  example  at  1  MHz  and  higher  we  don’t  expect  temperature- 
driven  variations  in  tissue  parameters  to  be  fast  enough  to  be  included  in  the  acoustic 
wave  equation  calculations.  The  simulations  to  be  shown  in  Chapter  5,  run  over 
the  time  of  several  seconds,  do  show  that  the  variation  in  background  parameters 
(especially  a)  is  very  important  to  temperature  predictions  in  the  focal  region. 

This  analysis  is  distinct  from  analysis  of  spatially-varying  media  as  a  debilitating 
factor  in  time  reversal  arrays.  Propagation  in  the  ocean  over  long  enough  ranges 
to  allow  the  channel  to  change  appreciably  over  the  time  scale  of  propagation  from 
source  to  array  is  possible.  Dowling  [22]  and  Dowling  and  Jackson  [21]  studied  this 
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effect,  ajid  conclude  that  while  the  time  reversal  process  is  robust  enough  to  allow 
for  good  retrofocusing  even  in  realistic  ocean  channel  conditions,  the  beamwidth  is 
increased.  Experimental  evidence  for  this  was  provided  by  Kuperman  et  al.  [52],  who 
suspended  a  400  Hz  vertical  line  array  in  the  Mediterraneaji  Sea,  and  were  able  to 
successfully  focus  onto  a  source  6.3  km  away  from  the  array. 
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Chapter  3 

FOCUSED  SOURCES  AND  ARRAYS 

In  this  chapter  basic  terminology  and  concepts  of  focusing  are  presented.  The 
application  of  focusing  to  time  reversal  arrays  and  therapeutic  ultrasound  devices  is 
left  until  later  specialized  chapters  on  these  topics.  The  means  of  controlling  the 
position  of  a  focal  spot  are  described  here,  and  some  simple  implementations  of  the 
wave  propagation  model  developed  in  Chapter  2  are  shown  and  compared  to  solutions 
by  other  methods  for  verification. 

3.1  Fields  in  Free  Space 

Our  goal  is  to  use  multiple  sources  to  construct  an  extended  focused  source  or  an 
array.  We  begin  with  a  description  of  the  acoustic  field  due  to  a  single  source  much 
smaller  than  a  wavelength,  commonly  known  as  a  point  source. 

3.1.1  The  Green’s  function  in  free  space 

A  periodically-varying  mass  in  unbounded  space  having  spherical  symmetry  has  a 
harmonic  complex  field.  The  Green’s  function  at  some  spatial  location  x  due  to  a 
source  at  x^,  denoted  by  ^(xlx^),  is  defined  as  the  solution  to  the  inhomogeneous 
Helmholtz  equation. 


(V^  -t-  P)  G(x|xs)  =  -47r  ^(x  -  X5), 


(3.1) 
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where  the  Dirac  delta  function  5(x  —  x^,)  is  an  impulsive  source  at  position  Xj.  The 
solution  to  this  equation  in  an  unbounded  medium  is  the  free-space  Green’s  function, 

pikR 

G(x|x,)  =  — ,  (3.2) 

where  the  displacement  R=  |a:  —  a^sj.  This  shows  that  a  field  due  to  a  point  source  in 
free-space  will  experience  a  1  /R  decay  without  any  medium  attenuation.  A  focused 
source  uses  the  combined  fields  from  many  small  (or  an  extended)  sources.  The 
sources  are  geometrically  arranged  so  that  a  stronger  field  results  at  locations  where 
the  arrivals  from  the  individual  component  sources  meet.  The  gain  obtained  by 
allowing  the  fields  of  many  small  sources  (or  the  integrated  extended  source)  to  add 
up  is  called  geometric  gain.  For  linear  acoustics,  this  gain  is  the  result  of  superposition, 
and  is  an  algebraic  sum  of  the  individual  fields  of  the  component  sources. 

3.2  Focusing  Using  Arrays 

This  section  describes  an  acoustic  source  consisting  of  more  than  one  active  element 
(an  array)  driven  in  some  deliberate  fashion  to  affect  a  localized  region  of  heightened 
average  acoustic  intensity  in  space,  time,  or  both,  called  a  focus.  The  position  of  the 
focus  in  space  and  time  is  determined  by  the  geometry  and  phase  relation  of  the  array 
elements,  as  well  as  propagation  path  parameters. 

For  focusing  in  the  linear  acoustics  case  using  the  superposition  principle  and 
multiple  transducer  elements,  a  wide  body  of  literature  exists  demonstrating  the  the¬ 
ory  and  application  of  linear  array  beamsteering  [19,  62,  69].  The  focusing  can  be 
performed  mechanically  by  placement  of  the  elements,  exploiting  their  geometry,  or 
electronically  by  applying  phasing  to  the  driving  signals,  [62].  Lithotripsy  and  ther¬ 
apeutic  focused  ultrasound  treatment  are  examples  of  the  use  of  focusing  to  form  an 
intense  acoustic  field  at  some  desired  location,  while  diagnostic  imaging  and  sonar 
are  examples  of  using  electronic  beamforming  to  locate  and  identify  objects  in  the 


acoustic  beam. 
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3.2.1  Parameters  we  cannot  control 

In  most  applications  the  propagation  path  is  given  by  the  nature  of  the  application. 
In  underwater  and  biological  applications,  the  path  can  be  complicated.  The  presence 
of  scatterers,  three-dimensional  inhomogeneous  medium,  as  well  as  absorption,  non¬ 
linearity,  and  dispersion,  make  it  difficult  to  model  the  propagation  in  these  media. 
We  will  see  how  phased  arrays  and  time  reversal  systems  can  be  used  to  overcome 
some  of  the  difficulties  associated  with  obtaining  a  focus  in  such  complicated  media. 

3.2.2  Parameters  we  can  control 

Typically,  we  can  specify  the  geometry  of  an  array  only  during  its  manufacture.  Since 
the  array  is  usually  an  expensive  component  of  the  hardware,  arrays  must  be  carefully 
designed  for  the  purpose  at  hand  for  best  efficiency.  Once  an  array  is  chosen,  the 
phasing  of  the  elements  is  the  parameter  most  often  used  to  create  and  position  the 
focus.  This  process  is  called  beamforming,  and  has  long  been  used  to  advantage  in 
underwater  engineering  applications.  Mechanical  manipulation  of  the  array  as  well 
as  electronic  beamsteering  can  be  used  to  move  the  main  lobe  of  a  focused  array. 
Diagnostic  ultrasound  equipment  generates  images  of  extended  regions  using  such 
techniques  to  construct  a  multi-dimensional  image  from  several  snapshots  (B-scans) 
in  the  receive  mode  along  a  direction. 

3.3  Linear  Arrays 

Perhaps  the  simplest  manifestation  of  multiple-element  sources  is  the  linear  array, 
consisting  of  multiple  elements  arranged  along  a  straight  line.  Note  that  a  linear 
array  refers  to  its  geometry,  and  should  be  distinguished  by  context  from  the  term 
used  to  describe  small-signal  acoustics.  Figure  3.1  shows  a  diagram  of  a  simple  linear 
array.  For  an  inter-element  separation,  d,  and  a  sound  speed,  c,  we  can  calculate 
some  relations  for  phasing  the  elements  of  a  linear  array  for  focusing  purposes  by 


Figure  3.1:  A  line  array  of  sources  phased  by  time  delays 


controlling  the  time  delay  for  the  element,  At„.  The  richness  of  the  device  can  be 
appreciated  by  considering  the  following  few  control  scenarios: 

•  Beamsteering:  making  the  time  delays,  At„,  increase  linearly,  At^+i  =  At„  +  ^. 
This  will  cause  the  main  wave  front  to  propagate  along  a  direction  other  than 
normal  to  the  face  of  the  transducer  array.  The  time  delays  to  affect  this  steering 
are, 

At„  =  n—sinO  +  T.  (3.3) 

c 


•  Focusing:  including  a  quadratic  variation  in  the  time  delays  to  obtain  focusing 
at  some  range,  F,  as  well  as  steering  angle,  9.  In  this  case  the  time  delays  are 
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given  by  Von  Ramm  and  Smith  [62]  as, 


Af„  = 


F 

c 


(3.4) 


•  Sidelobe  reduction:  including  array  shading,  whereby  various  amplitude-weighting 
windows  are  applied  to  the  array  aperture.  This  is  also  known  as  apodization. 
The  cross-sectional  intensity  profile  of  a  beam  may  be  controlled  using  such 
techniques.  This  can  be  used  to  reduce  the  normally  unwanted  sidelobes  or 
grating  lobes  of  a  beam,  where  acoustic  energy  is  heightened  at  periodic  angu¬ 
lar  directions  by  constructive  interference  from  periodic  array  sources. 

Note  that  an  overall  time  offset,  T,  is  added  to  the  delays  to  preclude  any  of  the  time 
delays  in  this  causal  system  from  being  negative. 

Other  details  such  as  individual  element  directivities  are  not  considered  in  this  dis¬ 
sertation,  but  will  affect  the  overall  array  field  pattern.  In  the  case  of  linear  acoustics, 
such  field  patterns  are  routinely  studied  and  calculated  and  measured  for  individual 
arrays  at  the  time  of  manufacture.  For  finite-amplitude  acoustics  it  is  generally  not 
possible  to  calculate  field  patterns  due  to  the  nonlinearity  and  possible  interaction 
between  the  fields  of  the  individual  elements,  except  via  approximations,  laboratory 
measurements,  and  computer  simulations. 


3.4  Arrays  with  Polar  Axial  Symmetry 

Figure  3.2  shows  some  geometries  of  focused  sources  with  axial  symmetry.  Such 
source  geometries  are  commonly  found  in  applied  and  research  settings.  In  (a)  a 
single  element  in  the  shape  of  a  bowl  achieves  its  focus  at  the  center  of  curvature. 
Single-element  spherical  bowls  are  common  in  hyperthermia  research  labs,  and  while 
they  provide  little  flexibility,  they  are  relatively  simple  to  use.  Another  simple  single¬ 
element  focused  source  is  shown  in  (b),  and  uses  an  acoustic  lens  with  a  different 
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sound  speed  than  the  surrounding  propagation  medium  to  produce  the  focusing.  This 
configuration  is  also  simple  to  use,  but  requires  knowledge  of  the  properties  of  the 
lens  and  the  propagation  medium  [54,  53].  In  addition,  problems  of  heating  can 
occur  in  the  lens  due  to  absorption.  Another  common  source  is  shown  in  (c),  and 


(a)  (b)  (c)  (d) 


Figure  3.2:  A  focused  source  can  be  obtained  by  manufacturing  a  single  element  as 
a  bowl  (a),  a  single  flat  piston  source  with  a  curved  acoustic  lens  (b),  an  array  of 
elements  arranged  on  a  bowl  (c),  or  a  flat  annular  array  (d). 


consists  of  multiple  small  elements  arranged  on  the  face  of  a  bowl-shaped  substrate. 
This  configuration  requires  sophisticated  driving  electronics  to  operate  correctly,  but 
provides  more  flexibility  in  forming  the  focus  than  a  single-element,  and  allows  more 
advanced  control  of  sidelobes  for  example.  Finally,  (d)  shows  a  flat  annular  array, 
which  allows  for  positioning  of  a  focus  along  the  axis  of  symmetry  of  the  array.  Grating 
lobes  in  annular  arrays  depend  on  geometry,  number  of  elements,  and  wavelength,  but 
can  be  reduced  by  increasing  the  bandwidth  of  the  signals  [23].  While  (a)  and  (b)  are 
only  capable  of  one  focus,  the  use  of  banks  of  array  elements  in  groups  allows  (c)  and 
(d)  to  form  multiple  foci  if  desired  [24]. 

Numerically,  we  consider  arrays  of  finite  aperture  as  consisting  of  many  electroa¬ 
coustic  transducer  elements  which  have  no  appreciable  spatial  extent  themselves,  and 
which  may  be  driven  independently  by  an  external  signal  generator.  This  suggests 
the  term  “sampled  aperture”  as  a  descriptor  of  an  array.  Another  common  modeling 
assumption  is  to  take  the  frequency  responses  of  all  the  elements  to  be  identical,  i.e., 
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the  elements  have  the  same  transfer  functions.  If  we  assume  that  the  total  resulting 
acoustic  field  is  the  linear  superposition  of  the  fields  from  each  of  the  small  sources 
we  can  compute  the  CW  field  quite  easily.  Similarly,  if  the  sources  are  pulsed,  then 
a  time-varying  version  of  the  field  is  obtained.  The  focusing  is  now  both  spatial  and 
temporal,  and  simple  time-of-flight  calculations  allow  us  to  compute  when  and  where 
the  maximal  intensity  region  will  occur.  Conversely,  the  focus  can  be  achieved  at 
the  desired  position  and  time  by  controlling  the  applied  pulses  with  time  delays  in  a 
predictable  manner. 

While  this  study  is  primarily  interested  in  the  time  domain  behavior  of  arrays, 
analogous  statements  can  be  made  for  the  frequency  domain  representation  of  the 
array  response.  The  frequency  response  of  the  array  is  the  superposition  of  the  fre¬ 
quency  responses  of  the  array  elements. 

An  array  may  be  used  in  a  transmit  (active)  mode  or  in  a  receive  (passive)  mode 
assuming  the  transduction  properties  of  the  elements  (both  electric  and  acoustic) 
have  such  symmetry.  This  allows  many  array  systems  to  be  used  for  transmission 
of  acoustic  energy  towards  a  given  position,  or  to  be  used  as  a  receiving  device  for 
recording  sound,  the  greatest  sensitivity  being  at  the  focus.  Also,  we  can  make  the 
distinction  between  focusing  at  some  range  and  beamsteering  along  some  steering 
angle,  however  we  will  refer  to  the  combined  effect  of  these  localizations  in  space  as 
well  as  in  time  simply  as  focusing. 

For  random  or  inhomogeneous  media,  the  Green’s  function  of  the  propagation 
medium  may  not  be  known  a  priori,  and  the  effects  of  medium  variability  may  not 
be  predicted  entirely  without  using  measuring  instruments  or  numerical  simulations. 

3.5  A  Numerical  Example:  1-D  Propagation 

We  now  present  a  simple  example  of  wave  propagation  using  the  model  developed  in 
this  chapter.  For  the  purposes  of  numerical  modeling,  the  nonlinear  term  is  expanded 
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into  two  terms,  and  the  absorption  coefficient  is  used.  The  wave  equation  is  discretized 
and  solved  in  the  following  form, 
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(3.5) 


Solving  an  equation  using  finite  differences  involves  discretizing  the  differentials  in 
the  partial  differential  equation  via  an  expansion.  The  continuous  variables  are  repre¬ 
sented  by  discrete  approximations,  given  in  a  matrix  representing  the  computational 
grid  and  stored  in  computer  memory.  The  finite-difference  time-domain  (FDTD) 
method  is  one  of  the  oldest  and  most  intuitive  methods  for  solving  differential  equa¬ 
tions  on  the  computer.  Details  of  the  numerics  are  left  for  the  Appendix. 

For  our  example,  a  one-dimensional  plane  wave  is  propagated  from  a  plane  wave 
source  in  a  homogeneous  medium  according  to  equation  (3.5),  in  Cartesian  coordi¬ 
nates, 
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(3.6) 


3.5.1  Solution  Using  the  FDTD  Method 

Solving  (3.6)  numerically  involves  discretizing  the  unknown  pressure  field,  p(x,  t)  onto 
the  spatial  and  temporal  grids  in  computer  memory.  The  discrete  representation  of 
the  pressure  is  p{xi,  tn),  or  simply  p",  with  integers  i  =  [1,2,...  ,Imax]  and  n  = 


,].  The  grid  locations  are  given  by 

Xi  —  Xq  d”  1)^37  5 

(3.7) 

in  —  ^0 

(3.8) 

where  5t  and  are  the  separation  between  adjacent  grid  locations  is  time  and  space. 
Note  that  the  present  discussion  is  only  valid  for  uniform  grid  spacing.  The  spatial 
extent  of  the  simulation,  or  the  computational  domain,  is  a:  =  [ro,  (/max  <^a;)],  and  the 
simulation  time  frame  is  t  =  [to,  (Amax  <^t)]- 


The  representation  of  the  differentials  is  derived  in  the  Appendix,  but  we  list  here 
the  relevant  expressions  for  the  terms  in  the  wave  equation,  accurate  to  0{5^,  5^): 

_ _ 

^5 _ 

^  =  ^'(3p?-4pf'+pr")+0(i?), 

^2 

^  =  4(Pr'^‘  -2p?  +  P?-*')  +  0(Sl),  (centered)  (3-9) 

_ ^4 _ 

=  ^  (2pf  -  +  4jor”^  - pr^)  (nght-sided) 

^3 

^  ^  (6p”  -  23pr^  +  34pr'  -  24pr"  +  -  P?~')+0(S^). 

We  assume  that  we  know  what  the  initial  and  all  previous  conditions  of  the  pressure 
field  are,  i.e.  all  variables  are  known  except  the  future  values  of  p(x,  t).  All  pressures 
prior  to  t  =  0  are  initialized  to  zero  in  our  simulations.  A  driving  pulse  is  applied  at 
source  locations  on  the  grid,  which  propagates  according  to  the  wave  equation.  The 
expressions  (3.9)  are  substituted  into  (3.6)  and  rearranged  to  solve  for  pressure  at  the 
next  time  step,  , 

-^1  +  92  ^3  +  92  ^P?  ^4  +  -  ^2,  (3.10) 


where 
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S.5.£  Boundary  Conditions 

We  impose  absorbing  boundary  conditions  (ABC)  at  the  edges  of  the  computational 
domain  to  prevent  or  to  minimize  reflections  from  the  edges  of  the  domain.  By 
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doing  so  the  simulations  approximate  the  behavior  of  waves  in  infinite  media.  Low 
order  ABC’s  are  relatively  easy  to  implement,  for  example  using  Kosloff  and  Kosloff’s 
technique  which  applies  a  progressively  more  absorbing  layer  near  the  edges  of  the 
domain  [50],  while  Mur’s  well-known  method  [60]  applies  a  radiation  condition 


dp'  1  dp' 
dx  Co  dt 


(3.12) 


The  ABC’s  and  their  implementation  will  be  left  for  the  Appendix. 

For  simulations  possessing  polar  cylindrical  symmetry,  only  one  half  of  the  two- 
dimensional  (r,  x)  domain  is  calculated,  or  r  =  [0,rTOax])  x  =  [xminiXmax]-  The 
boundary  condition  at  r  =  0  is 


dr 


(3.13) 


The  boundary  conditions  along  the  other  edges  use  radiation  condition  ABC’s. 
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Chapter  4 

TIME  REVERSAL  ARRAYS 

Several  factors  are  cited  for  reducing  the  quality  of  focusing  or  imaging  with  phased 
arrays  [62]: 

•  nonideal  response  of  transducer  arrays  and  limitations  in  delay  line  systems 

•  refraction  errors  due  to  medium  inhomogeneities 

•  target  ambiguities  due  to  phase  quantization  errors 

A  linear  time  reversal  array  (TRA)  or  phase-conjugate  array  can  correct  for  all  but 
the  last  deficiency  of  phased  array  systems  listed  above.  By  so  doing,  the  TRA  uses 
the  medium  as  a  matched  filter  to  automatically  generate  the  delays  at  each  element 
for  transmission  and  focusing  [20].  Time  reversal  arrays  are  alternately  called  time 
reversal  mirrors  (TRM)  owing  to  the  way  in  which  the  array  acts  as  a  temporal  mirror 
for  the  captured  signals. 

Time  reversal  mirrors  have  found  several  uses  in  medical  and  underwater  applica¬ 
tions  [72,  67,  52].  The  concept  of  time  reversal  is  an  extension  of  phase  conjugation 
theory,  which  is  known  to  hold  for  linear  fields  in  reciprocal  media  [48,  27].  The  time 
reversal  array  captures  not  only  the  initial  phase,  given  by  At„,  but  also  the  waveform 
by  sampling  signals  coming  from  the  target  through  a  propagation  medium  which  ide¬ 
ally  remains  stationary.  The  basic  operation  of  a  TRM  is  explained  in  a  number  of 
articles  by  Fink  [27,  28],  and  is  illustrated  in  Figure  4.1  as  a  three-step  process  consist¬ 
ing  of  target  illumination,  signal  collection,  and  time-reversed  retransmission.  Most 
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time  reversal  studies  consider  an  abbreviated  sequence  of  events,  starting  with  the 
array  receive  mode  (b),  followed  by  the  array  transmit  mode  (c). 


Figure  4.1:  Time  reversal  array  operation  after  Fink  [27].  In  (a)  the  target  is  illumi¬ 
nated  by  a  pulse  from  a  single  element,  in  (b)  the  receive  mode  where  the  illuminated 
target  acts  as  a  small  scatterer  emitting  a  spherical  pulse  while  the  array  is  recording, 
in  (c)  the  transmit  mode  the  array  retransmits  the  time-reversed  pulses  recorded  in 
the  receive  mode. 


4.1  Impulse  Response  Analysis 

Fink  [27]  has  described  the  behavior  of  linear  time  reversal  systems  in  terms  of  the 
impulse  response.  The  impulse  function  S{t)  can  be  used  to  describe  the  dynamics 
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of  the  reception  and  retransmission  of  acoustic  pulses.  The  impulse  response  analysis 
to  follow  is  after  Fink  [27].  If  an  array  element  Ei  causes  a  (non dimensional)  velocity 
excitation  v{t)  =  5{t)^  the  acoustic  velocity  potential  at  some  point  in  space,  Xo,  is 
described  by  the  diffraction  impulse  response  /j*(xo,t)  for  the  illuminating  (transmis¬ 
sion)  event  due  to  Ei.  A  related  impulse  response  /i^(xo,t)  can  be  defined  for  the 
receiving  event  due  to  an  impulse  source  at  Xq.  A  very  useful  aspect  of  the  time 
reversal  system  is  that  the  exact  impulse  responses  do  not  need  to  be  known  or  cal¬ 
culated.  What  is  necessary  is  that  the  transmit  and  the  receive  modes  have  identical 
impulse  responses,  whatever  these  may  be. 

This  analysis  is  dependent  on  acoustic  reciprocity.  For  a  Green’s  function,  G(xo,  to\x,  t), 
representing  the  field  at  x  at  time  t  due  to  an  impulse  excited  at  Xq  at  time  to,  the 
reciprocity  theorem  for  the  Green’s  function  states  that 

G(xo,  to|x,  t)  =  G(x,  to|xo,  t).  (4.1) 

We  can  write  the  diffraction  impulse  responses  in  terms  of  the  Green’s  functions  over 
the  transducer  surface,  S,, 


/i^(x,f)  = 

/  G(xo,fo|x,f)da:, 

JSi 

(4.2) 

hl{xo,t)  = 

f  G{x,to\xo,t)dx. 

(4.3) 

JSi 


Therefore,  we  need  not  distinguish  between  the  receive  and  the  transmit  mode  diffrac¬ 
tion  impulse  responses,  since  they  are  identical  by  (4.1).  Using  the  notation  of  Fink 
[27],  we  simply  denote  these  responses  by  /ij(xo,t). 

Another  set  of  transfer  functions  need  to  be  accounted  for.  That  is  the  transfer 
functions  describing  the  conversion  of  acoustic  energy  into  electrical  energy  during 
the  receive  mode,  h“®(t),  and  the  conversion  of  electrical  energy  into  acoustical  energy 
during  the  transmit  mode,  h®“(f).  These  transfer  functions  are  convolved  with  those 
of  the  propagation  during  the  time  reversal  process.  The  steps  of  the  time  reversal 
operation,  starting  with  emission  of  the  illumination  pulse  from  the  target  are; 
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1.  The  receive  mode:  For  a  target  located  at  Xo,  the  observed  (illumination)  signal 
is  the  convolution 

(g)hi(xo,t).  (4.4) 

The  emission  from  the  illuminated  target  is  sampled  by  the  TRM  over  a  time 

interval  T,  which  we  expect  to  contain  the  desired  signal  information.  The 
time-reversed  signal  is  thus 

hr{T-t)®hi{xo,T-t).  (4.5) 

2.  The  transmit  mode:  The  time-reversed  signals  must  now  be  retransmitted  by 
the  array,  requiring  convolution  through  /i?“(t)  and  propagation  to  the  target 
through  hi{xo,t), 

hT{T  -t)®  hi{xo,  T  -  t)  ®  hr(t)  (g)  hi{xo,  t).  (4.6) 

The  maximum  response  is  obtained  at  time  t  =  T,  and  if  =  hf^(^t),  which  is  a 

reasonable  assumption  for  high-quality  transducers.  The  reason  is  that  the  maximum 
output  of  a  linear  system  having  impulse  response  h{t)  is  achieved  when  the  input  is 
of  the  form  h{—t).  Such  a  response  is  calculated  from  the  convolution  h{t)  ®  h{—t), 
and  is  an  even  function  with  a  maximum  at  t  =  0. 

The  total  field  due  to  the  array  is  then  the  sum  over  the  elements, 

/i“®(T  -t)®  hi{xo,  T  -t)®  h|“(t)  ®  hi{xo,  t).  (4.7) 

i 

The  elements  achieve  constructive  interference  at  x  =  Xo  and  t  =  T,  resulting  in  a 
strong  focal  spot  field  akin  to  that  due  to  geometrical  focusing. 

4.2  Time  Reversal  Array  (TRA)  Debilitating  Factors 

One  goal  of  the  present  study  was  to  examine  the  feasibility  of  using  TRAs  in  a  shal¬ 
low  water  environment  to  focus  acoustic  energy  onto  waterborne  mines  for  possible 
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neutralization  by  acoustic  pressure.  Due  to  the  resolution  and  dynamic  range  lim¬ 
itations  of  current  systems,  uncertainties  in  amplitude,  time-domain  waveform,  and 
phase  will  result.  This  will  cause  a  general  degradation  of  the  performance  of  the 
phased  array  system.  An  effective  mine  countermeasure  (MCM)  system  will  certainly 
require  going  to  such  high  intensities  that  nonlinear  properties  of  the  transduction 
devices  and  the  propagation  medium  will  become  important.  The  effect  of  finite- 
amplitude  propagation  on  the  performance  of  a  TRA  and  the  effect  of  absorption 
in  the  propagation  medium  axe  also  studied  as  debilitating  effects  on  TRA  focusing. 
These  debilitating  effects  will  be  studied  in  the  remainder  of  this  chapter. 

4.3  TRA  Initial  Phase  Error  Effect 

Error  in  the  initial  time  delay  for  TRA  elements  is  the  first  debilitating  factor  we 
will  examine  in  this  chapter.  This  can  be  expected  to  be  one  of  the  most  important 
debilitating  factors  which  degrade  the  performance  of  a  TRA  system.  Since  the  TRA 
is  a  phased  array,  the  integrity  of  the  temporal  or  phase  information  is  crucial  to 
its  operation.  The  most  important  phase  data  is  that  which  determines  the  initial 
arrival  times  of  the  pulses  from  the  array  elements  at  the  location  of  a  target.  This 
feature  is  shared  with  ordinary  phased  arrays  or  time-delay  focusing,  which  relies  on 
time  of  flight  of  the  pulses  to  form  a  focus  by  superposition  at  the  target  location. 
Other  phase  effects  would  involve  the  phase  information  for  the  pulses  emitted  by  each 
TRA  channel  following  initial  activation,  this  would  determine  the  focusing  quality 
due  to  multipaths,  multiple  scattering,  and  phase  aberration  correction  from  medium 
refraction  index  inhomogeneities. 

4.S.I  Signal  Phase  Error 

Error  in  the  signal  initial  phase  was  modeled  by  introducing  random  stochastic  jitter 
into  the  time-domain  signals  that  the  array  elements  record  during  the  receive  mode. 
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The  jitter  can  result  from  limited-resolution  data  acquisition  electronics  in  the  receive 
and/or  transmit  modes.  The  jitter  is  calculated  for  each  array  element  individually 
as  a  time  delay  in  the  initial  phase  of  each  signal.  Because  the  present  study  seeks 
to  define  a  starting  point  for  evaluating  the  relative  effects  of  the  more  important 
parameters  described  that  reduce  the  effectiveness  of  the  time  reversal  array,  only 
the  simplest  cases  are  studied.  For  example,  one  could  have  introduced  jitter  into 
each  cycle  of  the  wavetrains,  or  into  each  time  step.  In  addition,  amplitude  jitter 
could  be  included.  The  simulations  in  this  study  only  introduce  jitter  into  the  initial 
phase  information  of  the  pulses  because  the  latter  effects  should  be  investigated  in 
the  category  of  pulse  shape  uncertainties,  which  is  not  addressed  here. 

The  jitter  is  given  in  terms  of  a  fraction  of  the  narrow-band  signal’s  period.  For 
a  time-domain  signal  at  array  element  k  of  the  form  we  introduce  a  time  delay, 
5tk,  so  that 

Pk{t)  Pk{t  -|-  Stk)-  (4-8) 

The  base  source  waveform  is  a  sinusoidal  envelope.  The  uncertainty  is  introduced  for 
each  of  the  elements  independently,  padding  the  leading  (jitter)  time,  5t^  with  zeros. 
The  jitter’s  duration  is  computed  randomly  for  the  element  from  the  narrow-band 
period,  tq,  the  maximum  error.  A,  as  a  fraction  of  2ir  of  the  base  wavelength  for  a 
run,  and  a  random  multiplier,  crjt, 

Stk  =  cTkAto.  (4.9) 

The  random  variable,  CTk,  can  range  from  zero  to  unity,  and  is  different  for  each 
element,  but  the  maximum  possible  jitter  for  any  array  element  is  A  for  a  given 
simulation.  Of  course,  the  jitter  can  be  defined  in  other  ways,  and  could  be  thought 
of  as  being  due  to  two  processes:  one  during  the  receive  mode,  and  the  other  during 
the  transmit  mode  of  the  array. 
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4.3, 2  TRA  Simulations  in  a  Shallow  Water  Channel 

A  vertical,  64-element,  equally-spaced  linear  array  with  an  aperture  of  25.6  m  is 
located  in  the  center  of  a  shallow  water  channel.  A  2  kHz  narrow-band  point  source 
is  located  51.2  m  away  from  the  array,  shown  in  Figure  4.2. 


Figure  4.2:  Diagram  of  the  layout  of  the  channel  and  time  reversal  system  simulated 
in  this  study.  The  block  diagram  of  the  time  reversal  system  is  after  Jackson  and 
Dowling  [48]. 


The  Solution  Method 

The  linear  inhomogeneous  acoustic  wave  equation  was  used  for  the  simulations  pre¬ 
sented  in  this  section. 


1  d'^p 


(4.10) 


\  p  J  pc^  dt^ 

with  the  primary  variable  being  the  acoustic  pressure,  p(r,f). 

The  wave  equation  is  solved  in  the  time  domain  using  a  two-dimensional  second- 
order  accurate  (in  space  and  time)  finite-difference  time-domain  (FDTD)  code.  The 
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calculations  are  carried  out  on  a  rectangular  grid  space  of  dimensions  1024  horizontal 
by  512  vertical  mesh  points.  Absorbing  boundary  conditions  were  used  at  the  extreme 
upper  and  lower  edges  of  the  computational  domain  to  simulate  an  extended  spatial 
region  for  visual  clarity,  although  this  is  not  necessary,  as  the  time  reversal  method 
is  valid  in  situations  where  multiple  paths  and  scattering  exist. 

The  Propagation  Medium 


SSP(x,z) 


x(m) 


Figure  4.3:  The  sound  speed  profile  used  in  the  simulations.  The  extreme  values  for 
sediment  and  air  are  truncated  to  better  illustrate  the  profile  in  water. 


The  wave  equation  (4.10)  is  solved  in  an  inhomogeneous  medium  modeled  as  air 
above  a  water  channel  approximately  50  m  deep  with  a  (fast)  fluid  sediment  below  it. 
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The  sound  speed  profile  (SSP)  is  fashioned  after  data  given  in  [25]  and  [1].  A  graphical 
representation  of  the  SSP  is  shown  in  Figure  4.3.  Although  the  basic  SSP  profile  is 
similar  to  that  used  in  many  studies,  it  serves  to  illustrate  the  physical  concepts 
only,  and  is  not  meant  to  be  an  accurate  oceanographic  representation  of  the  SSP 
of  any  actual  body  of  water.  The  properties  used  for  20°C  air  were  a  homogeneous 
sound  speed  of  343  m/s  and  a  density  of  1.2  kg/m^.  The  sound  speed  in  the  water 
channel  was  a  function  of  depth  below  the  surface,  with  inhomogeneities  added  to 
that  profile.  The  sediment  also  had  inhomogeneities  built  on  top  of  an  average  sound 
speed  of  1650  m/s,  and  an  average  density  of  1860  kg/m^.  Spectral  statistics  of  the 
channel  properties  were  not  considered  for  the  present  study.  The  density  field  was 
obtained  in  a  similar  manner  to  complement  the  SSP.  Inhomogeneities  in  the  water 
and  sediment  are  in  the  form  of  small  fluctuating  regions  of  excess  sound  speed  and 
density.  Figure  4.3  gives  an  idea  of  the  roughness  of  this  scale  and  its  magnitude. 
Furthermore,  a  flne  random  component  is  added  to  the  sound  speed  and  density  to 
give  some  fine  structure. 

The  shapes  of  the  air-water  and  water-sediment  interfaces  are  composed  of  combi¬ 
nations  of  sinusoids  with  small,  local,  random  fluctuations.  Again,  the  intention  is  to 
provide  simulations  in  a  non-uniform  medium  with  rough  interfaces  and  not  to  model 
any  oceanographic  spectra  at  this  time. 

4-3.S  Results  for  TRA  in  shallow  water  channel  simulations 

The  results  of  the  simulations  confirm  that  the  time  reversal  method  successfully 
focuses  linear  acoustic  waves  onto  a  target  in  an  lossless  inhomogeneous  medium  with 
multipath  effects  and  rough  boundaries.  Two  snapshots  are  shown  for  the  reference 
(ideal)  case  run.  Figure  4.4  shows  the  pressure  field  during  the  receive  mode  of  the 
array  in  the  top  panel,  and  the  instant  of  maximum  focus  onto  the  source  during  the 
transmit  mode  of  the  array  in  the  lower  panel. 
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Figure  4.4:  The  acoustic  pressure  field.  The  top  frame  shows  the  pressure  some  time 
after  leaving  the  source  during  the  array’s  receive  mode  operation.  The  lower  frame 
shows  the  pressure  after  the  array’s  transmit  mode,  when  the  maximum  pressure 
occurs  at  the  source  location. 


Description  of  TRA  Phase  Jitter 

Random  time-domain  jitter  is  introduced  in  the  form  of  zero  padding  leading  the 
initial  phase  signal  from  each  array  element.  Figure  4.5  shows  the  pressure  field  in 
dB  around  the  location  of  the  source  at  the  instant  of  maximum  focusing.  The  values 
are  referenced  to  the  maximum  pressure  (at  the  source’s  location  usually).  It  was 
found  that  the  location  of  the  focus  maximum  remained  near  the  original  location  of 
the  source.  The  reason  is  that  for  many-element  arrays  the  focal  shift  would  tend  to 
average  out  to  its  original  value  because  the  error  has  a  zero  mean.  On  the  other  hand, 
significant  degradation  in  focus  quality  was  observed  for  jitter  exceeding  about  one- 
tenth  of  a  wave  period.  As  expected,  the  initial  phases  of  the  waveforms  were  shown 
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to  be  very  important  to  the  focusing  ability  of  the  axray.  Simulations  were  conducted 
with  jitter  that  ranged  from  10%(2  tt)  up  to  a  full  27r  of  a  period.  The  -3  dB  points 
did  not  show  any  appreciable  spreading  from  case  to  case,  but  the  magnitude  of  the 
acoustic  pressure  for  the  cases  with  large  jitter  was  far  reduced  and  spread  over  a 
large  region  of  the  channel,  resulting  in  significant  focusing  degradation  for  jitters 
greater  than  10%  to  20%  of  one  period  (see  Figure  4.5).  This  result  was  published 
in  [41]  and  is  confirmed  for  line  arrays  in  general  by  Von  Ramm  and  Smith  [62],  who 
state  that  a  maximum  phase  error  equating  to  A/8  is  tolerable  for  arrays  with  at 
least  16  elements.  Wang  et  al.  [70]  also  find  that  for  a  two-dimensional  hyperthermia 
applicator  a  phase  error  equivalent  to  A/6  is  tolerable.  In  the  studies  cited  above,  the 
authors  find  that  an  acceptable  focus  was  maintained  for  random  phase  errors  up  to 
15  to  20%  of  a  wavelength. 

The  surface  and  bottom  reflections  can  act  as  virtual  sources  to  enlarge  the  ef¬ 
fective  aperture  of  the  time  reversal  array  [5].  Further,  the  inhomogeneities  and  the 
multipath  reflections  can  act  to  enlarge  the  effective  aperture  of  the  array.  The  latter 
phenomenon  has  been  recognized  by  Dowling  and  Jackson,  who  refer  to  it  as  “super- 
focusing”  in  [22].  The  localized  differences  in  the  index  of  refraction  act  as  distributed 
sources  which  are  present  throughout  the  medium.  The  focal  zone’s  full  width  at  half 
maximum  (FWHM)  spot  size  in  free-space  is  given  by  the  diffraction  limit  as 

w  =  1.2  \z/a,  (4-11) 

the  width  is  proportional  to  the  range,  .2,  and  inversely  proportional  to  the  aperture, 
a.  In  our  case  the  FWHM  would  be  approximately  1.8  m  in  free-space.  We  observe 
that  the  FWHM  points  in  the  channel  occur  at  ±  0.3  m  from  the  source  location, 
and  attribute  this  result  to  the  virtual  array  due  to  the  multipath  propagation  and 
the  superfocusing  described  above. 
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Figure  4.5:  Slices  through  the  source  position  at  the  time  of  maximum  focusing  for 
various  jitter  conditions.  The  dashed  lines  denote  the  vertical  slices  (parallel  to  the 
array),  while  the  solid  lines  denote  the  horizontal  slices  (perpendicular  to  the  array). 
Panels  axe  for  (a)  No  jitter  (reference  case),  (b)  Max.  jitter  10%(27r),  (c)  Max.  jitter 
20%(27r),  (d)  Max.  jitter  =  100%(27r). 


4.3.4  Phase  Jitter  Simulation  Conclusions 

The  premise  of  using  phase  conjugate  arrays  for  the  focusing  of  intense  acoustic  fields 
onto  a  remote  waterborne  target  has  implications  for  MCM  system  design.  The 
possibility  of  remote  neutralization  of  pressure-sensitive  mines  would  be  an  asset 
to  the  MCM  arsenals  that  exist  today  [31].  The  concept  has  been  demonstrated  in 
theory  and  in  the  laboratory  for  ultrasonic  frequencies  in  medical  applications.  Ocean 
field  experiments  in  shallow  water  have  been  conducted  recently  by  Kuperman  et  al. 
[52].  These  experiments  were  conducted  at  400  Hz,  and  would  not  encounter  the 
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difficulty  with  electronic  jitter  that  high  frequency  arrays  would  suffer.  However, 
the  results  obtained  so  far  from  experimental  [18]  and  theoretical  [22]  groups  show  a 
remarkable  robustness  when  using  the  time  reversal  technique  with  multiple  scattering 
and  reflection  in  random  media. 

We  used  linear  acoustics  to  model  sound  propagation  through  a  shallow  water 
channel  with  an  inhomogeneous  sound  speed  profile  as  well  as  a  rough  surface  and 
bottom.  The  results  for  cases  that  are  expected  to  degrade  the  focusing  ability  of  a 
time  reversal  array  by  altering  the  initial  phase  information  are  given.  Initial  phase 
timing  is  corrupted  by  some  jitter  in  the  time-domain  signals,  which  affect  the  relative 
phases  of  the  transmitted  array  element  waveforms. 

The  array’s  focusing  appears  to  hold  up  well  under  these  circumstances  for  jitter 
up  to  10%  to  20%  of  the  narrow-band  period.  Jitter  greater  than  about  20%  of  a 
period  results  in  significant  loss  of  focal  pressure.  These  results  are  confirmed  by 
other  studies  of  linear  phase  quantization  errors  [62,  70].  The  location  of  maximum 
pressure  remained  at  the  location  of  the  source  because  the  jitter  is  a  zero-mean 
random  variable,  indicating  that  initial  phase  of  the  returned  signals  is  more  important 
than  the  details  of  the  waveform  phase  shape.  This  is  encouraging,  since  the  data 
aquisition  of  a  broadband  time-reversal  signal,  and  the  translation  of  that  signal 
into  a  corresponding  high-intensity  array  pressure  output  is  unlikely  with  current 
technology.  This  is  especially  true  if  the  array  consists  of  elements  whose  bandwidth 
is  significantly  greater  than  PZT  transducers. 

Other  factors  not  presented  thus  far  which  are  detrimental  to  the  focusing  of  phase 
conjugate  arrays  are  the  nonlinear  behavior  of  the  medium  and  the  transducer  and 
electronics’  transfer  functions.  The  nonlinearity  of  the  propagation  process  plays  an 
important  role  in  conjunction  with  absorption  that  will  be  studied  in  Section  4.5,  and 
casts  serious  doubt  on  the  feasibility  of  using  time  reversal  systems  for  real  MCM 
applications. 
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4.4  The  Effect  of  Absorption  on  TRA’s 

In  this  section  the  linear  absorption  of  the  propagation  medium  will  be  examined 
as  a  debilitating  factor  on  the  performance  of  a  focusing  TRA.  The  motivation  for 
including  absorption  as  a  debilitating  factor  is  that  some  applications  for  TRA  focus¬ 
ing  may  be  carried  out  in  media  where,  unlike  water,  the  absorption  coefficient  may 
significantly  affect  the  propagating  waves. 

One  example  of  absorbing  media  is  biological  tissue.  This  section  presents  an 
example  of  a  TRA  in  a  coupling  water  bath  used  to  focus  acoustic  energy  onto  a 
target  inside  a  modeled  human  head.  The  transcranial  focusing  experiences  a  drastic 
effect  on  time  reversal  when  crossing  the  highly-absorbing  bone  layer  of  the  skull. 


4-4-i  Violation  of  Time  Reversal  Invariance  by  Absorption 


We  use  a  linear  version  of  the  model  wave  equation  (2.21),  derived  in  Chapter  2, 


vV- 


1  d^p 


-Vp  •  Vpo  + 


2a  d^p 


=  0, 


(4.12) 


Po  '  COO)*  afi 

where  the  absorption  term  is  written  in  terms  of  the  absorption  coefficient,  a,  related 
to  the  acoustic  diffusivity  5  by 


a  = 


2(? 


(4.13) 


If  we  consider  the  left  hand  side  (LHS)  of  (4.12)  as  a  linear  differential  operator, 
£,  acting  on  the  acoustic  pressure,  we  can  write  (4.12)  as 


Up  =  LHS(4.12)  =  0.  (4.14) 

We  can  test  whether  the  differential  operator,  £,  is  time-reversible  by  examining  its 
effect  on  the  forward  time  solution  and  its  backward  time  counterpart.  If  p  =  <;^(x,<) 
is  a  forward  time  solution  of  (4.12),  that  is,  </>(x,  f)  satisfies 


£<^(x,  t)  =  0, 


(4.15) 
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then  by  direct  calculation  we  can  test  whether  q{x,t)  =  —t)  is  a  solution  of  the 

wave  equation  as  well.  If  so,  then  (4.12)  is  time  reversal  invariant,  i.e.,  it  holds  under 
time  reversal. 

Since  for  t  =  —t,  the  partial  derivative  dr/dt  =  —  1,  upon  making  the  transfor¬ 
mation  t  T,  the  time  derivatives  of  p  are 

aXx,-i)  _  5>(x,t)  f  dry  _  ^  ■,,„^Xx,t) 

df-  ~  dr-  \dtj  ^  ’  dr^  '  ’ 


That  the  odd-ordered  partial  derivatives  with  respect  to  time  undergo  a  sign  change 
is  significant  in  our  discussion.  This  simple  fact  is  correctly  cited  for  holding  the  key 
to  the  time-reversibility  of  the  lossless  acoustic  wave  equation,  and  has  been  noted  by 
Fink  [27]. 

Now  we  substitute  for  p(x,t)  its  time-reversed  companion  solution,  q(x,  i),  into 
(4.12)  and  use  (4.16)  to  obtain 


c2  dt^  ) 


Po 


(x,-t) 


(4.17) 


from  which. 


Cq  —  jC(j>  — 


28d^<i> 
c4  dt^ ' 


(4.18) 


Thus  we  conclude  that  q(x,  t)  =  ^(x,  —t)  is  not  a  solution  of  (4.12)  in  the  presence  of 
absorption.  Physically  speaking,  this  says  that  the  linear  wave  equation  (4.12)  is  not 
time  reversible  due  to  the  presence  of  the  absorption  term. 


Example  of  TRA  in  an  Absorbing  Medium 

We  use  a  numerical  example  to  illustrate  the  detrimental  effects  absorption  has  on 
time  reversal  focusing  systems.  Figure  4.6  shows  the  scenario  presented  in  this  ex¬ 
ample,  which  involves  a  two-dimensional  simulation  of  the  acoustic  field  due  to  a  line 
array  operated  as  a  TRA  in  a  water  coupling  bath  to  insonate  a  scatterer  in  the  human 
brain.  A  scatterer  is  positioned  at  some  location  within  a  cross-section  of  a  human 
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Figure  4.6:  Diagram  of  the  propagation  model,  showing  the  position  of  the  linear 
array  to  the  left  of  the  head,  and  the  point  scattering  target  within  the  brain. 


head  containing  various  inhomogeneous  organs,  and  most  importantly,  a  bone  skull 
layer.  The  propagation  properties  for  the  organs  and  the  skull  bone  were  obtained 
from  the  literature  [77,  55,  30].  Ellipses  were  used  to  generate  the  geometrical  layout 
of  the  head  cross-section,  in  a  fashion  similar  to  that  used  to  generate  tomographic 
phantoms  [49,  55].  Twenty-three  ellipses  were  used  to  construct  the  2-D  phantom  in 
Figure  4.6.  A  line  array  of  64  elements  were  simulated  to  lie  in  the  water  about  1 
cm  from  the  head.  Simulations  of  TRA  focusing  with  and  without  absorption  were 
performed  and  the  results  on  the  focal  pressure  amplitude  were  compared. 
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A  second-order  FDTD  solution  for  the  acoustic  pressure  was  solved  on  the  2-D 
computational  domain  to  simulate  an  illuminated  point  scatterer  emitting  a  pulse 
while  the  array  collected  data  in  the  receive  mode.  This  step  was  followed  by  trans¬ 
mittal  of  unamplified  time-reversed  signals  from  the  array  elements,  which  will  focused 
onto  the  target.  Due  to  the  limited  aperture  provided  by  the  array,  not  all  the  il¬ 
lumination  pulse  energy  was  captured  by  the  array.  Hence,  the  time-reversed  focus 
amplitude  was  only  a  fraction  of  the  original  amplitude,  even  in  the  lossless  case. 
Figure  4.7  shows  the  time  traces  of  the  acoustic  pressure  as  measured  at  the  64  array 


Figure  4.7:  Time  traces  measured  by  each  of  the  64  TRA  elements  during  the  receive 
phase  of  the  time  reversal  operation. 


elements.  To  lowest  order,  the  initial  delay  time  was  due  to  time-of-flight  time  rep¬ 
resenting  the  distance  separating  the  array  element  and  the  illuminated  target  point. 
These  were  the  signals  which  will  be  time  reversed  by  the  TRA  during  the  transmit 
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phase  of  the  process. 

Figure  4.8  shows  the  acoustic  pressure  along  the  x  and  the  y  axes  taken  through 
the  target  location  at  the  instant  of  maximum  focusing.  The  waveform  emitted  by 


Figure  4.8:  Spatial  slices  of  the  pressure  field  along  the  x  and  the  y  axes  passing 
through  the  target  at  the  instant  of  maximum  focusing. 


the  target  and  the  waveform  returned  to  the  target  should  be  the  same,  except  scaled 
in  amplitude  and  time-reversed.  Figure  4.9  shows  the  original  pulse  sent  out  by  the 
scatterer  during  the  receive  mode  in  (a),  and  the  time  trace  at  the  target  location 
upon  time  reversal  focusing  from  the  transmit  phase  (b).  Note  that  some  of  the  high- 
frequency  features  have  been  lost  in  the  two-way  propagation.  This  is  because  the 
pulse  used  as  the  illumination  signal  was  obtained  from  laboratory  measurements  with 
broadband  instruments  of  a  shock  wave  generated  by  a  spark  lithotripter  source,  and 
the  absorption  in  the  fluid  acted  as  a  low-pass  filter.  For  the  case  where  thermoviscous 
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Figure  4.9:  Original  pressure  pulse  as  sent  out  by  the  scatterer  (a),  and  the  pulse  as 
measured  at  the  target  upon  time  reversal  (b)  for  the  case  with  absorption. 


absorption  was  included,  the  focal  amplitude  was  reduced  due  to  linear  absorption  due 
to  thermal  and  viscous  losses,  as  explained  in  Chapter  2,  and  time  reversal  invariance 
violation,  as  explained  in  the  previous  section.  Figure  4.10  shows  a  comparison  of 
peak  pressure  at  the  target  location  at  the  time  of  maximum  focusing.  The  spatial 
slices  through  the  target  location  at  the  instant  of  maximum  acoustic  focusing  pressure 
show  that  the  lossless  case  provided  a  much  higher  amplitude  at  the  focus.  In  our  case 
the  peak  pressure  at  the  focus  for  the  lossless  case  was  2.2%  of  the  original  pressure  of 
the  scatterer  pulse.  When  absorption  was  taken  into  account  the  pressure  at  the  focus 
dropped  to  1.3%  of  the  maximum  pulse  pressure.  This  reduction  in  focal  pressure  by 
about  one  half  was  due  to  the  two  effects  mentioned  in  the  previous  paragraph.  The 
reduction  in  focal  pressure  is  not  due  to  propagation  impediments  such  as  enhanced 
acoustic  impedance  mismatch  for  example.  The  same  sound  speeds  and  densities  are 
used  for  both  lossless  and  lossy  simulations.  The  severe  inhomogeneity  (the  skull 
layer)  was  responsible  for  the  fact  that  focal  pressure  was  so  low  compared  to  a  case 
where  no  absorption  would  occur.  For  this  reason  a  section  of  skull  bone  is  normally 
removed  surgically  from  the  array  aperture  prior  to  insonation  of  the  brain  in  animal 
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Figure  4.10:  Acoustic  pressure  at  the  focus  at  the  instant  of  maximum  focusing, 
showing  the  lossless  case  achieving  almost  twice  the  peak  pressure  as  the  lossy  Ccise. 


experiments.  Otherwise  the  large  absorption  in  the  skull  bone  would  cause  excessive 
heating  near  the  skull.  This  aspect  of  biomedical  acoustics  will  be  discussed  in  detail 
in  Chapter  5,  which  deals  with  hyperthermia  applications. 

Finally,  we  present  snapshots  of  the  acoustic  pressure  for  both  the  lossless  and 
the  lossy  simulations.  Figure  4.11  shows  one  snapshot  of  the  pressure  field  during  the 
receive  mode  of  the  TRA  operation,  where  (a)  is  the  lossless  case  and  (b)  is  the  lossy 
case.  Then  a  pair  of  frames  are  taken  from  the  transmit  stage  of  the  TRA  operation 
at  the  instant  of  maximum  focusing.  Frame  (c)  is  the  field  for  the  lossless  simulation, 
and  frame  (d)  is  the  field  at  the  same  time  for  the  lossy  simulation.  While  the  gray 
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Figure  4.11:  Snapshots  of  the  acoustic  pressure  during  the  receive  mode  (a)  and  (b) 
and  the  transmit  mode  at  focusing  (c)  and  (d).  The  frames  on  the  left  (a)  and  (c) 
are  for  the  lossless  simulation,  while  the  frames  on  the  right  (b)  and  (d)  are  for  the 
lossy  simulation. 


scale  was  fully  utilized  in  each  of  the  four  frames,  the  scales  to  the  side  show  that 
the  peak  pressure  for  the  lossy  case  was  about  one  half  of  that  of  the  lossless  case. 
Also  it  can  be  seen  from  the  transmit  mode  frames  and  from  Figure  4.10  that  the 
full-width-at-half-maximum  (FWHM)  for  the  lossy  case  is  increased.  In  this  example, 
the  FWHM  went  from  4  mm  for  the  lossless  case  to  9  mm  for  the  lossy  case. 

We  have  seen  how  absorption  can  degrade  the  focusing  ability  of  a  TRA.  It  might 
seem  that  one  solution  to  partially  overcome  this  debilitating  effect  is  to  amplify  the 
signal  at  the  array.  However,  if  the  amplification  is  such  that  the  signal  is  of  finite 
amplitude  we  encounter  further  debilitating  effects  due  to  nonlinearity,  which  will  be 


discussed  in  the  next  section. 
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4.5  The  Effect  of  Nonlinearity  on  TRA ’s 

Nonlinear  propagation  in  absorbing  sound  channels  is  proposed  as  a  debilitating  factor 
to  time  reversal  array  system  performance.  Since  we  have  seen  the  detrimental  effect 
of  linear  absorption  on  TRA’s  in  the  previous  section,  we  investigate  the  combined  role 
of  absorption  and  nonlinearity  in  this  section.  Here  we  ask  the  question;  what  extra 
role  if  any  does  nonlinear  propagation  of  finite-amplitude  sound  play  in  degrading  the 
time  reversal  invariance  of  a  TRA? 

4.5.1  Background  for  Nonlinear  TRA  Analysis 

We  now  look  at  the  effect  of  finite- amplitude  (nonlinear)  propagation  on  time  rever¬ 
sal  systems  operated  in  lossless  and  in  lossy  media.  A  finite-amplitude  sound  wave 
distorts  as  it  propagates  through  the  medium.  The  distortion  is  manifested  as  a  trans¬ 
fer  of  energy  from  the  originial  frequency  spectrum  into  higher  harmonics.  Because 
most  media  exhibit  absorption  that  increases  with  frequency,  one  might  expect  that 
nonlinear  distortion  would  lead  to  greater  absorption,  and  hence  worse  time  reversal. 
Muir  et  al.  [59]  conducted  an  experiment  where  a  finite-amplitude  sinusoidal  pulse 
was  propagated  through  water  and  then  reflected  off  of  a  pressure-release  surface. 
For  their  pulse  the  pressure-release  surface  caused  an  amplitude  multiplication  by 
— 1,  which  for  a  sinusoidal  signal  is  similar  to  introducing  a  180®  phase  shift.  They 
found  that  the  pulse  undistorted  after  reflection,  which  suggests  that  time  reversal 
invariance  may  hold  in  the  presence  of  nonlinear  distortion.  In  other  words,  while  it 
was  not  a  time  reversal  operation,  the  experiment  showed  that  nonlinear  distortion 
effects  could  be  undone. 

The  setup  for  the  problem  under  study  is  as  follows:  a  source  emits  a  plane  wave 
acoustic  pulse  which  propagates  to  a  time  reversing  array  of  two  elements,  positioned 
on  either  side  of  the  source.  The  acoustic  field  recorded  by  the  array  is  time  reversed, 
and  then  retransmitted.  This  simple  setup  is  suitable  for  studying  the  effects  of 
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nonlinear  distortion  of  an  acoustic  field.  Of  particular  interest  is  the  debilitating 
effect  that  nonlinearity  may  have  on  a  time  reversal  system’s  ability  to  form  an  intense 
focus  near  the  original  source.  This  is  of  practical  interest  to  those  wishing  to  use 
time  reversal  systems  to  achieve  high-intensity  acoustic  fields  at  a  target,  such  as  in 
the  destruction  of  kidney  stones  by  lithotripsy  [67],  the  ablation  of  tumors  in  tissue 
[66],  or  the  remote  neutralization  of  mines  at  sea  [41]. 

Two  scenarios  were  studied:  first,  a  finite-amplitude  pulse  was  emitted  from  a 
source  and  propagated  outward  to  the  time  reversal  array  (the  receive  stage  of  the 
time  reversal  operation).  The  signal  was  then  time  reversed  and  retransmitted  from 
the  array  elements  (the  transmit  stage  of  the  time  reversal  operation).  We  label  this 
scenario  the  nonlinear-nonlinear  case  because  the  traveling  pulse  had  finite  ampli¬ 
tude  during  both  the  receive  and  the  transmit  stages  of  time  reversal,  and  hence, 
undergoes  nonlinear  effects  during  both  stages  of  the  process.  The  second  scenario 
considered  was  where  the  source  emitted  a  low-amplitude  pulse,  which  propagated 
without  appreciable  nonlinear  steepening  to  the  array  during  the  receive  mode,  fol¬ 
lowed  by  amplification  of  the  time-reversed  signal  at  the  array  during  the  transmit 
mode.  The  amplification  for  the  transmit  mode  propagation  led  to  nonlinear  dis¬ 
tortion  of  the  pulse.  This  case  is  referred  to  as  the  linear-nonlinear  scenario  in  our 
study. 


4.5.2  The  Nonlinear  Absorbing  Wave  Equation 

The  model  equation  used  in  this  study  for  propagation  of  finite-amplitude  acoustic 
pressure  pulses,  p(x,t),  in  a  thermoviscous  medium  was  derived  in  Chapter  2, 

1  d^\ 


Vpo  ^  s  d^p  /?  d'^p^  _ 
p - : —  Vp  -t-  — H - —  U. 


Po 


Cq  dfi  poCq 


(4.19) 


The  last  term  on  the  left  hand  side  of  (4.19)  is  the  nonlinear  term,  discussed  in 
Chapter  2. 

The  most  salient  features  of  pulse  propagation  in  an  absorbing  nonlinear  medium 
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Figure  4.12;  Comparison  of  the  pulse  used  in  the  present  study  as  emitted  from  the 
source  (solid  line)  and  as  observed  some  time  later  (dashed  line).  The  pulse  steepens 
due  to  nonlinearity  and  decays  due  to  absorption. 


are  the  amplitude  reduction  due  to  absorption  and  the  waveform  distortion,  or  shock¬ 
ing,  due  to  the  nonlinearity.  An  example  of  the  pulse  used  in  this  study  before  and 
after  propagating  some  distance  in  an  absorbing  nonlinear  medium  is  shown  in  Fig. 
4.12. 

4.5.3  The  Conditions  for  Time  Reversal  Invariance 

Following  the  method  used  in  the  section  4.4,  we  now  consider  the  left  hand  side  (LHS) 
of  (4.19)  as  a  nonlinear  differential  operator,  V,  acting  on  the  acoustic  pressure,  and 


write  (4.19)  as 


Vp  =  LHS(4.19)  =  0. 

If  /)  =  ^(x,t)  is  a  forward  time  solution  of  (4.19),  that  is,  <^(x,  t)  satisfies 

V(f){x,  t)  =  0, 


(4.20) 

(4.21) 


then  we  can  test  whether  g(x,  t)  =  ^(x,  -t)  is  a  solution  of  the  wave  equation.  If 
so,  then  (4.19)  is  time  reversal  invariant,  i.e.,  it  holds  under  time  reversal.  We  can 
expect  the  absorption  term  to  lead  to  time  reversal  invariance  violation,  as  was  found 
in  Section  4.4,  but  we  will  also  uncover  effects  due  to  nonlinearity. 

Substituting  for  p(x,t)  its  time-reversed  companion  solution,  q{x,t),  into  (4.19) 


we  obtain 


c2  p 


C4  I  dt^ 


pd^  df^ 


(4.22) 


from  which. 


Vq  =  T>(j) 


25d^(j> 

c4  dt^  ■ 


(4.23) 


Thus  we  conclude  for  the  nonlinear  absorbing  wave  equation  that  q(x,  t)  =  <^(x,  -t) 
is  not  a  solution  of  (4.19)  in  the  presence  of  absorption,  as  for  the  linear  absorbing 


wave  equation. 

Mathematically,  the  result  indicates  that  the  nonlinear  term  itself  is  not  respon¬ 
sible  for  time  reversal  violation.  However,  as  stated  earlier,  real  nonlinear  wave  prop¬ 
agation  is  always  accompanied  by  absorption  in  the  medium.  The  result  is  an  ir¬ 
reversibility  of  nonlinear  acoustics.  Therefore,  for  finite-amplitude  propagation  we 
would  expect  a  degradation  of  the  performance  of  a  time  reversal  system.  The  above 
discussion  applies  to  our  first  test  scenario  for  this  study,  the  nonlinear-nonlinear  case. 

Now  we  modify  the  analysis  for  the  second  scenario  under  study,  the  linear- 
nonlinear  case.  If  the  acoustic  pressure  received  is  not  only  time  reversed,  but  am¬ 
plified,  then  the  correct  test  for  time  reversal  invariance  checks  V  for  the  case  where 
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q{'x.,t)  =  a0(x,  —t),  where  a  is  a  constant  real  amplification  factor.  In  this  case  a<p  is 
substituted  for  (j>  in  (4.22), 


Vq 


c2  dt^  ) 


(3  d^{a<f>y 
pc*  df^ 


a 


V(t>- 


’pc*  dt^ 


(4.24) 

(x,-t) 

(4.25) 


Note  that  the  propagation  is  strictly  “nonlinear”  in  both  travel  directions,  but  the 
nonlinear  effects  are  negligible  during  the  first  (receive)  leg  of  the  travel  due  to  low 
amplitude.  The  amplitude  is  amplified  by  using  a  >  1  for  the  transmit  mode  of  the 
operation. 

The  result  above  shows  that  when  the  signal  undergoes  amplitude  amplification, 
time  reversal  invariance  in  a  nonlinear  medium  no  longer  holds,  even  if  the  fluid  is 
lossless.  For  time  reversal  invariance  to  hold  for  nontrivial  solutions  this  case  requires 


that 

1.  a  =  1  or  /?  =  0,  and 

2.  S  =  0. 


The  physical  explanation  is  that  aside  from  pure  scaling  of  the  solution,  the  change 
in  signal  amplitude  from  receive  to  transmit  modes  implies  that  different  nonlinear 
distortion  occurs  during  the  two  stages  of  the  time  reversal  process,  which  will  inval¬ 
idate  the  time  reversal  invariance.  While  not  strictly  a  violation  of  reciprocity,  this 
is  analogous  to  doing  time  reversal  in  a  medium  which  was  not  steady  for  the  dura¬ 
tion  of  the  two  stages  of  the  time  reversal  process.  The  receive  and  transmit  mode 
propagation  occur  in  what  is  essentially  different  filters  in  the  channel. 

4.5.4  Numerical  Study  of  Finite- Amplitude  TRA’s 

Numerical  simulations  were  used  to  investigate  the  behavior  of  time  reversal  systems 
while  controlling  the  absorption  coefficient  and  the  coefficient  of  nonlinearity  in  the 
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wave  equation.  The  peak  pressure  amplitude  in  the  receive  mode  (at  the  source)  and 
in  the  transmit  mode  (at  the  array)  was  independently  controlled  also  to  allow  for 
one-way  or  two-way  finite-amplitude  propagation,  as  desired. 

4.5.5  Description  of  the  Nonlinear  TRA  Simulations 

The  wave  equation  (4.19)  was  solved  numerically  for  the  acoustic  pressure,  p,  using  a 
finite-difference  time-domain  (FDTD)  code.  The  finite  differencing  was  fourth  order 
accurate  in  space  and  second  order  accurate  in  time.  The  simulations  were  performed 
for  a  homogeneous  thermoviscous  fluid.  Numerical  considerations  such  as  stability 
and  numerical  dispersion  limited  the  parameter  space  which  could  be  covered  using 
the  explicit  FDTD  method  used.  However,  good  resolution  of  the  acoustic  pressure 
field  was  possible  for  the  range  of  parameters  used  in  this  study. 


X  (cm) 


Figure  4.13:  Geometry  of  the  plane  wave  time  reversal  simulations  showing  a  temporal 
snapshot  in  the  array  receive  mode.  The  source  is  at  location  B,  while  the  two  array 
elements  are  at  locations  A  and  C. 
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4.5.6  Run  Geometry  and  Parameters 

The  simulations  were  carried  out  in  a  one-dimensional  coordinate  system  with  x 
denoting  the  spatial  coordinate.  The  source  from  which  the  source  pulse  was  emitted 
is  located  at  x  =  0.  Fig.  4.13  shows  the  location  of  the  source  and  the  elements 
with  relation  to  each  other.  The  source  element  (B)  was  driven  with  a  sinusoidal 
burst,  shown  in  Fig.  4.12,  having  a  duration  of  6  fxs,  a  center  frequency  f  —  I 
MHz,  and  a  Gaussian  envelope.  The  driving  pressure  waveform  applied  at  the  source 
had  a  maximum  pressure  amplitude  of  1  MPa  for  the  simulations  requiring  nonlinear 
receive  mode  propagation,  and  1  Pa  for  those  requiring  linear  propagation  in  the 
receive  mode  (where  the  computational  domain  was  much  shorter  than  the  shock 
formation  distance). 

The  calculations  presented  in  this  study  were  for  a  fluid  having  a  small  signal 
sound  speed  c  =  1500  m/s,  and  density  p  =  10^  kg/m^,  but  otherwise  was  not 
meant  to  represent  water.  The  values  chosen  for  the  absorption  and  the  nonlinearity 
coefficients,  as  well  as  the  frequency  and  the  peak  pressure  of  the  pulses  cover  a 
parameter  space  well-suited  for  the  numerical  method  used  in  the  study.  Effects  such 
as  absorption  and  nonlinearity  are  cumulative,  and  will  manifest  themselves  in  the 
waveforms  over  some  distance  and  time.  In  this  study,  we  attempted  to  demonstrate 
the  effects  under  study  over  a  span  of  50  wavelengths  using  full  wave  simulations 
suitable  for  detailed  study  of  the  pulses  in  space  and  time.  The  values  of  a  and  /?  are 
not  those  of  any  specific  fluids,  and  are  chosen  to  have  an  effect  on  the  propagation 
over  the  distance  used  in  the  simulations.  We  could  have  alternatively  used  the  value 
of  peak  acoustic  pressure  to  vary  the  shock  formation  distance,  for  example. 

.  In  order  to  simplify  and  generalize  the  results  some  useful  quantities  which  were 
discussed  in  Chapter  2  are  repeated  here  for  convenience.  The  acoustic  Mach  number 

€  =  \  (4.26) 

pc^ 

where  po  is  a  characteristic  acoustic  pressure  of  the  pulse.  The  shock  formation 
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distance,  Xghock,  is  the  distance  a  sinusoidal  plane  wave  travels  before  developing  a 
shocked  profile,  and  is  defined  by 


1 


(4.27) 


where  k  is  the  wave  number  given  by  2;r// c.  Note  the  dependence  of  e,  and  hence 
^shock,  on  the  peak  pressure  amplitude,  po-  Finally,  the  Gol’berg  number, 

r  =  4.  ('••28) 

aX 


is  used,  where  A  =  27r/A;. 


4.5.7  Results  for  Finite-Amplitude  TRA’s 

We  have  seen  in  Section  4.5.3  that  the  absorbing  nonlinear  wave  equation  treated 
as  a  differential  operator  acting  on  the  acoustic  pressure  should  be  time-invariant 
when  the  absorption  term  is  negligible  for  the  nonlinear-nonlinear  scenario.  For  the 
linear-nonlinear  scenario,  we  also  saw  that  the  absorption  has  to  be  negligible,  but 
in  addition,  that  either  the  amplification  at  the  array  must  be  negligible  (a  =  1)  or 
the  nonlinearity  coefficient  must  be  negligible  (/?  =  0).  In  this  section  we  show  the 
relative  importance  of  the  nonlinearity  and  the  absorption  in  time  reversal  using  the 
two  scenarios: 

1.  Nonlinear-Nonlinear  time  reversal:  In  this  scenario  the  source  is  driven  with 
a  high  peak  amplitude  (1  MPa),  resulting  in  nonlinear  steepening  occurring 
during  the  array  receive  mode.  The  detected  signals  at  the  array  elements  are 
time  reversed  without  modification  to  the  amplitude  of  the  signals  at  the  array 
elements.  The  retransmitted  pulses  leave  the  array  with  a  significant  fraction  of 
the  original  source  amplitude,  and  nonlinear  effects  occur  during  the  transmit 
stage  of  the  time  reversal  as  well.  The  values  of  the  absorption  coefficient,  a,  and 
the  nonlinearity  coefficient,  /?,  were  controlled,  and  were  the  same  for  a  given 
run  for  both  the  receive  and  the  transmit  stages  of  the  time  reversal  operation. 
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2.  Linear-Nonlinear  time  reversal:  A  low-amplitude  (1  Pa)  pulse  was  emitted  from 
a  plcine  source.  The  pulse  was  then  captured  by  two  array  elements  positioned 
on  either  side  of  the  source.  The  received  signals  were  time  reversed  and  retrans¬ 
mitted  from  the  array  elements  (also  plane  sources),  but  the  amplitude  of  the 
retransmitted  signal  was  amplified  at  the  array  by  scaling  it  up  to  1  MPa  peak 
pressure.  Again,  various  values  of  a  and  /?  were  used  to  investigate  the  effects  of 
time  reversal  in  nonlinear  absorbing  media.  For  this  scenario  the  low  amplitude 
of  the  receive  mode  propagatipn  ensured  that  the  outgoing  pulse  propagated 
without  undergoing  nonlinear  distortion,  as  the  shock  formation  distance  was 
much  larger  than  the  propagation  distance.  The  array  transmit  stage,  on  the 
other  hand,  caused  the  pulse  to  undergo  significant  distortion  due  to  the  large 
pressure  amplitude  provided  by  the  array. 


4-5.8  Nonlinear-Nonlinear  Time  Reversal 

This  case  may  apply  if  either  a  large  target  is  interrogated  or  illuminated  by  a  high- 
amplitude  pulse,  or  if  a  source  emits  a  high-amplitude  pulse  to  be  time  reversed.  Here 
we  expect  waveform  steepening  during  both  stages  of  the  time  reversal  process,  as 
the  signal  will  have  appreciable  amplitude  during  both  the  receive  and  the  transmit 
modes.  Note  that  capturing  and  retransmitting  a  shocked  waveform  experimentally 
requires  broadband  array  elements  and  electronics  to  obtain  good  fidelity. 

A  reference  simulation  was  run  in  which  both  the  absorption  and  the  nonlinearity 
were  suppressed  by  setting  a  =  0  and  /?  =  0.  This  run  gave  data  for  the  “ideal” 
time  reversal  case,  as  would  be  expected  given  our  eaxlier  analysis.  The  peak  positive 
pressure  obtained  at  the  source/target  location  upon  time  reversal  is  denoted  by 
Ptargef  waveforms  emitted  by  the  source  were  fully  recovered  in 

shape  and  amplitude  back  at  the  source  position  following  the  time  reversal  operation. 
Note  that  Ptarget  nonlinear-nonlinear  scenario  was  not  the  same  as  Ptarget 
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Figure  4.14:  Nonlinear-nonlinear  time  reversal:  Greatest  positive  half  cycle  of  the 
pulse  at  the  target,  normalized  to  the  ideal  case  (a  =  0,/?  =  0),  for  different  values  of 
|5,  with  a  =  10  Np/m.  Note  that  some  numerical  error  is  noticeable  in  the  waveform 
for  the  case  /3  =  40. 


the  linear-nonlinear  scenario,  but  within  each  scenario,  pressures  were  normalized  to 
the  peak  target  pressure  for  the  ideal  case  at  hand,  a  =  0  and  /?  =  0. 

Fig.  4.14  shows  the  time  trace  of  the  largest  positive  half  cycle  at  the  target  upon 
time  reversal,  normalized  to  Ptarget-  For  the  case  of  «  =  10  Np/m,  we  found  that  the 
best  retrofocusing,  with  /?  =  0,  resulted  in  only  47%  of  Ptarget-  When  /?  was  increased 
such  that  the  source-to-array  separation,  a:src-array,  was  1  shock  formation  distance,  the 
peak  target  pressures  achieved  were  about  37%  of  Ptarget-  The  peak  target  pressure 
decreased  monotonically  to  about  23%  of  ptarget  for  a  ^src-array  of  about  1.8  or 


65 


Figure  4.15:  Nonlinear-nonlinear  time  reversal:  Greatest  positive  half  cycle  of  the 
pulse  at  the  target  for  different  values  of  /?,  with  a  =  20  Np/m. 


a  of  40. 

For  a  higher  absorption  value,  a  =  20  Np/m,  the  pressure  at  the  target  was  even 
lower,  see  Fig.  4.15.  Here  the  peak  positive  target  pressure  iox  (i  =  0  was  about 
24%  of  Ptarget  absorptiou,  and  became  less  as  /?  was  increased  through 

50,  where  peak  pressure  was  only  13%  of  Pta^get  ^  source-array  separation  of  about 
2.3  ^Jshock* 

Fig.  4.16  shows  how  peak  pressure  at  the  target  varied  as  a  function  of  nonlinearity 
for  a  =  10  and  20  Np/m.  Nonlinearity  is  measured  in  terms  of  the  number  of  shock 


( 
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Figure  4.16:  Nonlinear-nonlinear  time  reversal:  Comparison  of  peak  pressure  at  the 
target  for  different  values  of  /?,  with  a  =  10, 20  Np/m.  While  source-to-array  distance, 
ajsrc-array  is  Constant,  increasing  /3  reduces  the  shock  formation  distance  aishock- 


formation  distances  separating  the  array  and  the  target.  This  can  be  related  to  by 

=  507r/?e.  (4.29) 

^shock 

It  was  noticed  that  the  rate  of  evolution  of  such  mature  pulses  slowed  down,  as  the 
creation  of  harmonics  by  nonlinearity  was  brought  into  balance  with  their  absorp¬ 
tion.  Also,  once  significantly  attenuated,  the  waveforms  did  not  experience  as  much 
nonlinear  distortion,  owing  to  their  smaller  pressure  amplitude. 

As  a  measure  of  the  effect  of  nonlinearity,  when  compared  to  the  case  /?  =  0,  we 
noted  a  relative  decrease  in  peak  target  pressure  of  about  22%  when  the  source-target 
separation  was  1  shock  formation  length  for  a  =  10  Np/m,  and  about  17%  decrease 
for  a  =  20  Np/m. 
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Figure  4.17:  Linear-nonlinear  time  reversal:  Greatest  positive  half  cycle  of  the  pulse 
at  the  target  for  dilferent  values  of  /?,  with  a  =  10  Np/m. 


4-5.9  Linear-Nonlinear  Time  Reversal 

This  scenario  is  the  most  common  in  acoustic  phase  conjugation  experiments,  where 
the  source  emits  (or  target  is  interrogated  with)  low-amplitude  pulses,  but  in  the 
transmit  mode  the  retransmitted  time  reversed  pulse  is  highly  amplified.  In  these 
simulations  the  source  emitted  a  1  Pa  acoustic  pulse,  which  was  captured  by  the 
array  elements  and  retransmitted  at  1  MPa,  an  amplification  factor  of  a  «  10®,  or 
120  dB. 

Fig.  4.17  shows  results  from  the  case  where  a  =  10  Np/m,  and  (S  was  varied 
between  0  and  40,  corresponding  to  P  between  0  and  7.5.  It  was  observed  that 
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Figure  4.18:  Linear-nonlinear  time  reversal:  Greatest  positive  half  cycle  of  the  pulse 
^  at  the  target  for  different  values  of  l3,  with  a  =  20  Np/m. 


in  all  of  the  runs  having  absorption  the  peak  pressure  at  the  target  following  time 
reversal  was  lower  than  Ptarget-  This  is  to  be  expected  from  our  proof  in  section  4.5.3. 
In  addition,  a  monotonic  decrease  in  retrofocusing  peak  pressure  was  observed  with 
increasing  nonlinearity  parameter.  Above  a  /?  of  40  the  waveforms  were  too  shocked  to 
be  simulated  using  the  current  method.  The  best  result  for  this  absorption  value  was 
obtained  for  /3  =  0,  in  which  case  the  peak  pressure  at  the  target  was  69%  of  Ptargef 
As  /?  was  increased,  equivalent  to  increasing  the  source-target  separation  up  to  about 
3.3  shock  formation  distances,  the  peak  target  pressure  was  reduced  to  about  51% 
of  Ptargef  while  the  absorption  was  clearly  responsible  for  time  reversal  focusing 
degradation,  the  presence  of  nonlinearity  and  the  amplification  at  the  array  further 
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reduced  peak  target  pressure  by  a  further  25%. 

More  dramatic  effects  were  seen  when  the  absorption  was  increased  to  20  Np/m. 
Fig.  4.18  shows  the  peak  pressure  half  cycle  at  the  target  for  the  higher  absorption. 
Here  ^  was  varied  from  0  to  70  before  numerical  errors  developed  in  the  highly  shocked 
waveforms.  As  before,  the  best  time  reversal  refocusing  was  obtained  when  (3  was 
zero.  In  this  case  the  peak  positive  pressure  at  the  target  was  50%  of  Pjargef  P 
wcis  increased  to  70  (a  F  of  6.5),  the  peak  attainable  pressure  dropped  to  less  than 
30%  of  Ptargef  The  effect  of  both  a  and  (3  on  the  peak  positive  achievable  pressure  at 
the  target  normalized  to  pj^arget  is  shown  in  Fig.  4.19  for  the  linear-nonlinear  scenario. 

As  a  measure  of  the  effect  of  nonlineaxity  in  the  linear-nonlinear  scenario,  when 
compared  to  the  case  =  0,  we  noted  a  decrease  in  peak  target  pressure  of  about  8% 
when  the  source-target  separation  was  1  shock  formation  distance  for  a  =  10  Np/m 
and  about  ’6%  for  a  =  20  Np/m.  This  was  about  one  third  the  relative  decrement 
noted  in  the  nonlinear-nonlinear  scenario. 

We  further  note  that  in  comparison  to  Figs.  4.14  and  4.15  where  approximately 
sinusoidal  waveforms  were  recovered.  Figs.  4.17  and  4.18  show  highly  distorted  pulses 
recovered  at  the  target.  This  is  because  in  the  linear-nonlinear  scenario  there  is  no 
nonlinear  distortion  from  the  receive  mode  to  be  undone  by  the  distortion  in  the 
transmit  mode,  recall  Muir’s  aforementioned  experiment  [59].  This  could  pose  a 
problem  in  applications  where  pulse  shape  reconstruction  is  a  measure  of  time  reversal 
system  performance. 

4.6  Conclusions 

This  chapter  described  time  reversal  arrays,  and  studied  the  effects  of  various  debili¬ 
tating  factors  on  the  ability  of  time  reversal  systems  to  form  intense  pressures  at  the 
focus  of  the  TRA.  Initial  phase  uncertainty  was  studied  in  a  simulated  shallow  water 
channel  containing  inhomogeneities  and  multipath  reflections.  It  was  found  that  the 
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Figure  4.19:  Linear-nonlinear  time  reversal:  Comparison  of  peak  pressure  at  the 
target  for  different  values  of  /?,  with  a  —  10, 20  Np/m.  While  source-to-array  distance, 
Xsrc-array  is  Constant,  increasing  /3  reduces  the  shock  formation  distance  ajshock- 


maximum  allowed  delay  errors  in  the  initial  signal  phasing  was  equivalent  to  about 
1/6  of  a  cycle  for  a  periodic  signal.  Most  applications  of  time  reversal  systems  in 
acoustics  exploit  their  ability  to  backpropagate  acoustic  energy  to  a  scattering  target 
or  a  source.  As  such,  it  is  important  to  know  the  limits  of  usefulness  for  such  retrodi- 
rective  focusing,  or  the  potential  for  achieving  maximal  intensity  at  the  focal  spot. 
One  issue  known  to  plague  time  reversal  systems  is  absorption  in  the  propagation 
medium.  Amplitude  compensation  has  been  used  with  some  success  to  correct  for 
amplitude  distortion  [66],  but  such  techniques  have  not  been  explored  for  nonlinear 
propagation,  in  which  case  the  amplitude  compensation  at  the  array  will  have  detri¬ 
mental  effects  of  its  own.  This  study  assumed  that  the  effectiveness  of  a  time  reversal 
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system  can  be  measured  solely  by  the  amplitude  of  the  pressures  at  the  target.  How¬ 
ever,  in  some  cases  the  waveform  shape  may  also  contribute  to  the  fidelity  of  the  time 
reversal.  Under  nonlinear  propagation  conditions  this  fidelity  will  be  compromised 
by  the  tendency  for  the  positive  portions  of  the  waveform  to  travel  at  greater  sound 
speed  than  the  negative  portions  of  the  waveform. 

It  was  shown  that  time  reversal  focusing  can  still  occur  in  absorbing  media,  since 
the  phasing  of  the  signals  is  not  affected,  however,  the  absorption  leads  to  a  violation 
of  time  reversal  invariance,  and  hence  the  effectiveness  of  the  TRA  is  compromised  if 
the  goal  is  achieving  a  high-intensity  focus.  While  the  spatial  extent  of  the  focal  spot 
is  unaffected  by  the  inclusion  of  absorption,  the  reduction  in  peak  amplitude  causes 
an  increase  in  the  FWHM  of  the  focus. 

Since  most  media  exhibit  a  frequency  power  law  of  absorption,  it  is  important 
to  know  what  effect  the  nonlinear  harmonic  generation  will  have  on  time  reversal 
systems  if  they  are  to  be  used  at  finite  amplitude.  Another  issue  raised  by  this 
study  is  how  amplification  of  the  signals  received  by  the  array  elements  will  com¬ 
promise  the  time  reversal  process.  We  showed  that  for  the  nonlinear-nonlinear  case 
the  nonlinearity  in  the  wave  equation  only  affects  the  time  reversal  invariance  in  that 
it  introduces  extra  absorption  due  to  the  generation  of  higher  frequency  content  in 
an  acoustic  waveform.  This  indirect  effect  appears  to  be  nearly  as  important  as  the 
medium’s  absorption  coefficient  in  cases  where  well-formed  shock  waves  appear,  or 
where  propagation  over  several  shock  formation  lengths  occurs.  Caution  is  in  order 
if  one  contemplates  increasing  the  gain  on  the  time  reversal  array  element  amplifiers 
in  an  attempt  to  overcome  the  focusing  loss,  as  in  the  linear-nonlinear  case.  Since 
the  shock  formation  distance  is  directly  related  to  the  peak  source  pressure,  shock 
formation  will  occur  sooner,  leading  to  a  reduction  in  the  performance  of  the  system 
for  a  given  source-array  separation.  We  conclude  that  altering  the  amplitude  of  the 
time  reversed  signals  at  the  array  will  lead  to  degradation  of  the  focus  in  a  nonlinear 
propagation  situation,  even  if  the  medium  is  lossless. 
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It  was  shown  that  for  propagation  over  several  shock  formation  distances  the  ef¬ 
fects  of  nonlinearity-enhanced  absorption  are  almost  as  important  as  the  effects  of 
the  absorption  itself.  The  extra  absorption  of  shocked  waveforms  will  significantly 
degrade  the  ability  of  a  time  reversal  system  to  deposit  a  high-intensity  acoustic  pres¬ 
sure  field  onto  a  target.  Additionally,  amplitude  alteration  at  the  array  elements  in 
a  nonlinear  medium  can  have  detrimental  effects  on  the  time  reversal  array’s  perfor¬ 
mance.  These  results  have  implications  for  future  uses  of  high-intensity  acoustic  time 
reversal  systems  such  as  mine  neutralization,  hyperthermia,  and  lithotripsy  using  time 
reversal  arrays. 
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Chapter  5 

TISSUE  HYPERTHERMIA  USING  FOCUSED  SOURCES 

5.1  Introduction 

There  exists  ample  motivation  to  explore  the  possibility  of  performing  noninvasive,  or 
bloodless,  treatment  of  deep-seated  tumors  in  the  human  body.  Traditional  surgery  to 
reach  and  remove  tumors  requires  cutting  a  path  to  the  affected  region,  controlling  the 
damage  to  the  intervening  organs,  controlling  the  resultant  bleeding,  and  maintaining 
infection-free  conditions  during  and  after  the  surgery.  The  concept  of  using  focused 
ultrasound  to  destroy  the  offending  tumor  tissue  by  heating  it  until  the  affected  cells 
die  due  to  the  thermal  dose  was  postulated  almost  half  a  century  ago  [29].  It  is  only 
recently  that  the  technology  to  build  and  control  and  image  the  focused  ultrasound 
surgery  has  become  possible.  These  efforts  have  intensified  in  the  past  ten  years,  as 
successful  experiments  in  animals  have  been  conducted,  and  a  small  but  aggressive 
industry  now  seeks  to  test  and  license  such  devices  for  treatment  in  humans. 

General  descriptions  of  focused  ultrasound  surgery  (FUS),  also  known  as  high- 
intensity  focused  ultrasound  (HIFU),  are  given  in  Hynynen  [47],  ter  Haar  [65],  and 
Sanghvi  and  Hawes  [63].  Research  on  therapeutic  ultrasound  is  interdisciplinary,  and 
is  of  interest  to  scientists  working  in  medical  applications  as  well  as  clinicians  with 
some  grasp  of  the  physics  and  engineering  of  ultrasound.  This  study  maintains  the 
perspective  of  physical  acoustics  and  does  not  emphasize  the  pathology  of  cell  injury. 
The  mechanism  for  the  destruction  of  the  afflicted  tissue  sites  is  known  as  necrosis, 
associated  with  the  pathologic  change  following  progressive  degradation  due  to  lethal 
injury  of  living  cells.  Cell  structure  and  membranes  are  known  to  suffer  as  a  result  of 
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thermal  overdose  from  the  deposition  of  acoustically-driven  heat  in  the  tissue.  Other 
mechanical  damage,  notably  as  a  result  of  cavitation  are  known  to  occur  at  high 
intensities  [46,  71],  but  are  not  investigated  here. 

Because  of  the  complexity  of  propagation  in  biological  tissue,  analytical  and  com¬ 
putational  solutions  to  pressure  field  and  heating  calculations  can  only  be  obtained 
using  simplifications  and  assumptions.  For  example,  when  describing  medical  ultra¬ 
sound  devices  and  their  effect  it  is  common  to  find  the  assumption  that  the  propaga¬ 
tion  medium  is  homogeneous  [9, 15],  linear  [36,  56,  58,  39],  or  both  [14,  24,  74,  32,  35]. 
Sophisticated  numerical  studies  have  been  conducted  recently  whereby  measured  tis¬ 
sue  properties  are  used  in  the  calculation  of  acoustic  fields  [57,  38,  40]. 

Despite  the  recent  advances  in  computational  power,  realistic  three-dimensional 
simulations  remain  slow,  costly,  and  conspicuously  absent  from  the  literature.  One 
reason  for  this  other  than  the  computing  shortage  is  that  3-D  measurements  in  tissue 
would  require  unconventional  measurement  techniques  such  as  MRI  or  tomography, 
and  would  be  much  more  difficult  to  obtain.  Some  important  problems  in  wave  prop¬ 
agation  can  be  studied  using  finite-difference  methods  in  2-D.  This  chapter  presents 
results  from  simulations  of  therapeutic  ultrasound  in  tissue-like  media.  The  questions 
addressed  are:  what  are  the  effects  of  nonlinearity  and  absorption  on  the  propagation 
and  heating  of  a  focal  region?  and  what  is  the  effect  of  medium  inhomogeneity,  and 
the  effect  of  temperature-dependent  tissue  parameters  on  the  heating  behavior  near 
the  focus  of  therapeutic  ultrasound  devices?  Single  pulse  and  continuous  wave  (CW) 
simulations  were  carried  out.  The  modeling  of  the  acoustic  propagation  and  focused 
sources  has  already  been  given  in  Chapter  2.  In  this  chapter  a  common  model  for 
the  thermal  behavior  of  tissues,  known  as  the  bioheat  equation,  will  be  explained, 
and  means  of  coupling  the  acoustic  and  the  thermal  aspects  of  the  problem  will  be 
presented. 

The  simulations  address  the  following  phenomena  and  their  effect  on  transient 
and  steady-state  acoustic  pressure  and  temperature  fields:  First,  the  nonlinear  effects 
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resulting  from  propagation  from  pulsed  devices  which  can  produce  finite-amplitude 
effects.  Second,  the  effect  of  inhomogeneities  in  two  dimensions,  as  deduced  from 
simulations  in  data  from  slices  of  human  tissue  obtained  by  experimental  measure¬ 
ments,  reported  by  other  researchers.  Third,  the  effect  of  temperature-dependent, 
time-varying  tissue  background  parameters.  It  is  found  that  inhomogeneity  of  soft 
tissue  can  have  the  effect  of  displacing  and  breaking  up  the  site  of  sound  and  heat 
deposition,  but  will  have  little  effect  on  the  overall  thermal  dose.  Finite-amplitude 
effects  are  more  important  in  calculating  the  size  and  temperature  of  the  hot  spot. 
Temperature-dependent  tissue  characteristics  appear  to  be  quite  important  in  predic¬ 
tion  of  thermal  effects  of  focused  ultrasound  from  CW  devices. 

5,2  Heating  Model:  The  Bioheat  Equation 


As  a  consequence  of  the  thermal  and  viscous  properties  of  fluids,  energy  loss  results 
when  an  acoustic  disturbance  passes  through  the  fluid.  The  acoustic  wave  deposits 
the  absorbed  energy  as  heat  [61,  14,  74].  For  the  case  of  tissue,  a  linear  bioheat 
equation  commonly  used  to  describe  the  thermal  effects  [16]  is 


dT  Ktiss  ■^2rp  ^bCb  ^rp  rp  \  ,  Q  I'l 

T  -  -pr-{T  -  r„)  4-  -^p-  (h-1) 

Ot  pOtiss  P^tiss  Pt^tiss 

where  T  is  the  difference  between  the  tissue  temperature  and  the  ambient  (arterial) 
temperature  (37°C'),  Kuss  is  the  thermal  conductivity  of  the  tissue  (0.6  W/m-K), 
Ctiss  is  the  heat  capacity  (3700  J/kg-K),  Wb  is  the  perfusion  or  cooling  by  blood  flow 
(0.5  kg/m^-s),  Cb  is  the  heat  capacity  of  the  blood  (3800  J/kg-K),  Ta  is  the  ambient 
arterial  temperature,  and  Q  is  the  heat  deposition  source  term,  described  in  Pierce 


|61], 
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(5.2) 


We  assume  that  a  is  the  loss  due  to  absorption,  and  scattering  loss  is  negligible  in 
the  attenuation  of  sound.  In  a  thermoviscous  fluid,  the  absorption  a  is  related  to  5 
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and  uj  =  2nf  by  [42] 


(5.3) 


5.2.1  FDTD  Solution  of  the  Bioheat  Equation 

In  a  manner  analogous  to  the  solution  of  the  wave  equation,  we  solve  (5.1)  using 
the  FDTD  method.  The  partial  derivatives  are  discretized  to  second  order,  and 
the  temperature  is  advanced  in  time  from  one  time  step  to  the  next.  Define  the 
variables  on  the  interior  grid  cells  having  spatial  and  temporal  discretization  6^  and 
5t  respectively. 
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The  partial  derivatives  are  written  to  second  order 

dt 


a^r 


AtVt, 
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dy 

and  the  source  term,  Q'^-  can  be  calculated  from 


Qh  =  A,’ 


(5.14) 


(5.15) 


The  bioheat  equation  then  can  be  used  to  solve  for  the  temperature  at  the  next  time 
step,  n  +  1, 
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5.3  The  Effect  of  Nonlinearity:  Single  Pulse  Simulations 


Theoretical  and  experimental  studies  of  nonlinear  propagation  and  focusing  in  tissue 
and  tissue-like  materials  have  shown  that  increased  heating  will  result  from  steepened 
or  shocked  waveforms  in  tissue.  There  has  been  increasing  interest  in  the  use  of  high 
intensity  ultrasound  in  tissue.  At  high  intensities  finite-amplitude  effects  can  lead 
to  the  production  of  nonlinearly  generated  harmonics.  These  harmonics  have  been 
exploited  in  recent  years  to  improve  imaging  capabilites  in  diagnostic  ultrasound 
machines  -  a  technique  commonly  referred  to  as  harmonic  imaging.  In  addition, 
focused  ultrasound  surgery  (FUS)  is  a  promising  technique  that  uses  high-intensity 
ultrasound  to  destroy  tissue  in  a  confined  region  which  will  avoid  the  necessity  of 
traditional  invasive  methods. 

The  explanation  for  increased  heating  from  nonlinear  ultrasonic  applicators  is 
explained  by  Wu  and  Du  [75]  and  Bacon  and  Carstensen  [3],  and  can  briefly  be  sum¬ 
marized  as  follows:  Most  materials  exhibit  a  frequency  power  law  for  the  absorption  of 
sound,  i.e.,  the  conversion  of  acoustic  energy  into  heat  energy.  For  single-frequency 
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waves,  the  small  signal  absorption  coefficient,  a,  depends  on  the  frequency  of  the 
acoustic  waves,  /,  and  a  real  exponent,  u, 

a~f\  (5.17) 

where  u  is  between  1  and  2,  and  is  about  1.1  for  most  soft  tissues  [2].  The  nonlinear 
generation  of  higher-frequency  harmonics  implies  excess  absorption  will  occur  for 
steepened  or  shocked  waveforms. 

While  bioeffects  have  been  widely  studied  in  the  laboratory  [65] ,  it  is  difficult  to 
obtain  detailed  spatial  measurements  without  disrupting  the  acoustic  field.  It  has 
been  known  for  some  time  that  high  intensity  ultrasound,  and  the  accompanying 
finite-amplitude  effects  result  in  different  bioeffects  than  linear  propagation  would 
predict  [3,  10].  In  a  thermoviscous  fluid,  the  absorption  is  proportional  to  the  square 
of  the  frequency,  so  as  higher  frequency  content  is  generated  during  wave  steepening, 
it  has  been  shown  that  increased  heating  will  result. 

5.3.1  Results  of  the  Nonlinear  Simulations 

Some  simulation  results  are  shown  in  this  section  to  illustrate  the  output  of  the 
numerical  code.  The  example  outputs  are  for  a  focused  bowl  source  array  of  64  simple 
sources  having  azimuthal  symmetry.  The  source  was  driven  by  a  1  MHz  sinusoidal 
burst  of  6  cycles  modulated  by  a  Gaussian  envelope  in  time.  The  geometric  focus  of 
the  array  was  situated  3  cm  from  the  array  face.  The  source  aperture  was  4  cm;  the 
computational  domain  spanned  an  area  of  5.12  cm  x  5.12  cm.  The  grid  spacing  was 
Ax  —  Ar  =  0.1  mm  and  At  =  10  ns. 

The  propagation  medium  was  modeled  as  a  tissue-like  material,  having  properties 
similar  to  those  reported  in  the  literature  for  soft  tissue.^®  The  parameters  used  in 
the  acoustic  problem  were;  cq  =  1600  m/s,  po  =  1100  kg/m^,  a  =  4.5  Np/m,  and 
(3  =  5.5.  For  the  bioheat  equation  the  following  parameters  were  used:  =  0.6 

W/(m-K),  Ct  =  Cb  =  3800  J/(kg-K)  and  Wt  =  0.5  kg/(m^-s).  Figure  1  shows  some 


79 


T(8(is) 


-2  0  2 


T(16|iO 


T(24tis) 


T(32ns) 

0 
1 
2 

3 

4 

5 

-2  0  2 


Figure  5.1:  Snapshots  of  the  pressure  (top  row),  and  temperature  (bottom  row)  for 
the  1  MPa  nonlinear  simulation.  The  axes  are  labeled  in  cm. 


snapshots  of  the  pressure  and  temperature  fields  resulting  from  nonlinear  propagation 
in  the  medium. 

Figure  5.3.1shows  axial  slices  of  the  temperature  fields  at  32  fis  for  linear  and 
nonlinear  simulations  at  source  pressures  ranging  from  1  MPa  to  10  MPa.  The  linear 
simulations  were  achieved  by  setting  (3  to  zero.  Temperature  elevation  increased  with 
source  pressure.  For  nonlinear  simulations  excess  heating  increased  dramatically  for 
source  pressures  above  5  MPa.  For  the  1  MPa  nonlinear  simulation  the  the  heating 
was  only  2%  above  the  linear  predictions  and  at  2  MPa  the  excess  heating  was  4%. 
However,  at  5  MPa  excess  heating  was  27%  and  at  10  MPa  it  was  80%.  High  source 
levels  were  required  to  observe  excess  heating  because  the  entire  propagation  path  was 
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(a)  (b) 


(a)  Axial  slices  of  the  temperature  at  32  yus  for  various  source 
pressures  with  and  without  nonlinearity,  (b)  Peak  temper¬ 
ature  elevation  for  several  source  conditions,  showing  the 
increasing  effect  of  nonlinear  distortion. 

modeled  as  tissue.  Tissue  has  a  relatively  large  absorption  and  prevents  steep  shocks, 
with  significant  harmonic  content,  from  forming.  For  propagation  through  water, 
which  is  much  less  absorbing,  steep  shocks  occur  much  more  readily  and  consequently 
excess  heating  becomes  important  at  lower  source  levels. 

5.3.2  Conclusions  from  the  Nonlinear  Effect  Simulations 

Accurate  knowledge  of  the  behavior  of  ultrasonic  beams  in  tissue  allows  for  better 
prediction  of  bioeffects  of  ultrasound,  and  improved  treatment  and  device  design.  The 
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FDTD  method  for  simulation  of  transient  finite-amplitude  acoustic  fields  has  been 
described  with  an  application  in  medical  ultrasound  used  as  an  example.  Heating 
in  tissue,  modeled  as  a  thermoviscous  fluid,  was  obtained  for  a  short  acoustic  pulse 
from  a  focused  bowl  source.  Peak  temperature  rises  from  a  1  MHz  ultrasonic  pulse 
propagating  nonlinearly  through  a  tissue-like  material  was  observed  to  depend  on  the 
degree  of  nonlinear  distortion  and  ranged  from  2%  to  80%  excess  temperature  rise  for 
1  MPa  to  10  MPa  source  pressure  conditions. 

5.4  The  Effect  of  Inhomogeneity:  A  2-D  Study 

In  this  section  we  look  at  the  role  played  by  tissue  inhomogeneity  in  the  context 
of  nonlinear  therapeutic  ultrasound  propagation.  A  1  MHz  pulse  of  ultrasound  is 
propagated  using  FDTD  simulations  in  2-D  Cartesian  coordinates.  The  propagation 
medium  was  constructed  from  2-D  measurements  of  tissue  propagation  delay  times 
obtained  from  researchers  at  the  University  of  Rochester,  described  by  Hinkleman, 
et  g/.,  [43,  44].  The  data  were  adapted  for  our  use  in  the  simulations  by  scaling  the 
magnitude  of  the  inhomogeneity  in  sound  speed,  density,  absorption  coefficient,  and 
nonlinearity  coefl&cient  to  the  desired  contrast.  In  the  original  study,  human  breast 
and  abdominal  wall  tissue  slices  were  preserved,  and  delay  times  were  measured  at 
fine  intervals  across  the  samples,  generating  a  2-D  map  of  the  inhomogeneities  in  the 
samples.  One  such  sample  is  used  in  this  study.  It  can  be  expected  that  local  effects 
such  as  diffraction  and  hot  spot  distortion  occur  in  inhomogeneous  media.  We  look 
at  this  effect  in  the  context  of  therapeutic  ultrasound  as  a  function  of  the  strength  of 
the  inhomogeneity. 

5.4-1  Description  of  the  Propagation  Medium 

The  inhomogeneity  contrast  (magnitude  of  the  variations  in  background  propagation 
parameters)  is  controlled  for  the  study,  and  varies  from  ±0%  (homogeneous  tissue) 
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Figure  5.2:  The  sound  speed  used  for  the  ±10%  inhomogeneity  contrast  case,  having 
a  base  value  of  1600  m/s  in  tissue  half  plane,  and  1500  m/s  in  the  water  half  plane. 
The  position  of  the  curved  source  and  geometric  focus  are  also  shown.  Note  that  the 
source  is  slightly  offset  to  the  left  to  avoid  the  occurrence  of  left-right  symmetries 
about  the  center. 


to  ±20%  contrast.  The  following  assumptions  are  made  in  generating  the  data  files 
for  the  background  properties: 

1.  Data  was  only  provided  for  delay  time,  or  sound  speed.  We  assume  the  inho¬ 
mogeneous  contours  for  sound  speed  delineated  regions  of  different  tissue  type 
and  composition.  Hence  it  is  assumed  that  all  tissue  properties  vary  along  the 
same  contour  lines  in  space. 

2.  The  percentage  of  contrast  in  each  of  the  background  properties  is  the  same. 
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i.e.  a  ±10%  variation  in  sound  speed  implies  a  ±10%  variation  in  density,  etc. 

3.  There  exists  a  qualitative  correlation  between  the  sign  of  the  variation  in  the 
properties  compaxed  to  the  time  delays.  The  assumption  used  here  is  that 
sound  speed,  density,  and  absorption  coefficient  have  the  same  tendency  to 
increase  or  decrease  together,  while  the  nonlinearity  coefficient  tends  to  violate 
this  direction,  i.e.  a  ±10%  variation  in  sound  speed  corresponds  to  a  ±10% 
variation  in  nonlinearity  coefficient.  This  last  assumption  is  borne  by  qualitative 
examination  of  tissue  measurement  data  given  in  the  literature  [33,  34,  77]. 

It  is  known  that  realistic  3-D  inhomogeneities  will  have  a  quantitatively  different 
effect  than  the  present  2-D  inhomogeneities,  but  3-D  data  is  not  available  for  such 
a  study  at  this  time.  We  expect  3-D  diffraction  to  have  a  more  severe  effect  on  the 
pressure  and  temperature  fields  than  2-D  diffraction,  and  so  the  results  reported  here 
are  expected  to  be  conservative  in  their  estimation  of  real  body  effects. 

The  inhomogeneous  initial  condition  data  file  for  a  medium  parameter,  Xoi^iV), 
which  is  a  function  of  2-D  space  is  generated  according  to  the  following  formula: 

T 

Xo{i,j)  =  Xoo  ±  Xoo  (5.18) 

where  xoo  is  a  reference  background  value  to  be  perturbed  according  to  the  template 
file,  Z,  which  contains  the  inhomogeneous  data  measured  and  normalized  to  lie  be¬ 
tween  [—1,  1],  and  r  is  the  contrast  (percent)  desired  for  the  run.  Only  one  half  of 
the  computational  domain  is  filled  with  this  inhomogeneous  tissue-like  medium.  The 
other  half  is  assumed  to  be  homogeneous  water  containing  the  focused  source.  An 
example  of  an  inhomogeneous  data  file  used  in  these  simulations  in  given  in  Figure 
5.2,  showing  a  ±10%  contrast  in  sound  speed.  The  measured  data  were  only  given  on 
a  32x128  grid  due  to  the  original  sample  size.  For  the  purposes  of  the  simulations, 
the  measured  data  was  extended  over  the  entire  tissue  portion  of  the  domain  by  tiling 
and  reflecting  it  to  fill  a  256x256  grid  region.  This  increased  the  apparent  size  of  the 
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Figure  5.3:  Acoustic  pressure  and  temperature  above  background  (37°C)  at  16/xs, 
24/us,  and  32/US.  The  water  half  plane  appears  unheated  due  to  the  very  low  absorption 
coefficient  in  water  compared  to  tissue.  Each  figure  is  scaled  independently  to  take 
advantage  of  the  full  gray-scale  color  range.  Axes  labeled  in  cm. 


tissue  sample,  while  minimizing  any  effects  of  repetition  on  the  sound  field,  e.g.  giving 
artificial  discontinuities,  which  could  cause  reflections.  Since  the  inhomogeneous  data 
is  symmetrical  about  the  x=1.28  cm  position  to  the  left  and  to  the  right,  the  source 
was  slightly  offset  to  the  left  so  that  the  left  and  the  right  portions  of  the  beam  do 
not  pass  through  identical  media. 
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Figure  5.4:  Acoustic  pressure  and  temperature  above  background  (37°C)  at  16)us, 
24^s,  and  32//s  for  the  ±10%  inhomogeneity  contrast  case.  Axes  labeled  in  cm. 

5.4-2  Results  From  the  Inhomogeneous  Medium  Study 

The  reference  (no  inhomogeneity,  zero  contrast)  run  was  one  having  two  half  planes, 
each  of  which  was  homogeneous.  The  water  half  plane  uses  the  following  values  for 
all  the  simulations  in  this  study: 

c  =  1500  m/s, 
p  =  1000  kg/m^, 
a  =  2.88  X  10"“*  Np/m, 

/3  =  3.5. 


(5.19) 
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The  mean  values  for  the  tissue  parameters,  about  which  the  inhomogeneous  files  were 
built  were: 

c  =  1600  m/s, 
p  =  1100  kg/m^, 

^  (5.20) 

a  =  4.5  Np/m, 

(3  =  5.5. 

These  parameters  are  similar  to  those  published  in  a  number  of  sources  over  the  last 
two  decades  [33,  34,  77].  Note  the  4  orders  of  magnitude  disparity  in  the  absorption 
coefficient  between  water  and  tissue.  For  this  reason,  the  temperature  profile  shows 
measurable  temperature  rise  only  in  the  tissue  half  plane,  while  the  water  half  plane 
heating  is  negligible. 

Snapshots  from  the  pressure  and  corresponding  temperature  calculations  for  the 
reference  run  are  shown  in  Figure  5.3.  Only  the  water-tissue  interface  interferes  with 
the  beam.  The  snapshots  are  taken  at  16,  24,  and  32  microseconds  from  launch  at 
the  source.  Peak  focusing  occurred  at  about  24  ps  as  the  pulse  reached  the  geometric 
focus.  The  thermal  hot  spot  continues  to  heat  up  even  after  24  ps  because  of  the 
slow  time  to  dissipate  heat  by  conduction  and  perfusion.  Upon  introducing  inhomo¬ 
geneities  within  the  tissue,  the  acoustic  focus  is  broken  up  and  distorted  at  various 
scales.  This  is  seen  in  Figure  5.4  along  with  its  effect  on  the  temperature  distribution 
for  the  ±10%  inhomogeneity  contrast  case.  Increasing  the  inhomogeneity  contrast  to 
±20%  results  in  the  fields  shown  in  Figure  5.5.  Slices  along  the  x  and  y  axes  are 
taken  from  the  acoustic  pressure  fields  for  the  contrast  cases  ±0%,  ±10%,  ±20%  at 
24/us  and  are  shown  in  Figure  5.6.  The  acoustic  field  is  scattered  and  distorted  by 
the  inhomogeneities  in  the  case  for  nonzero  contrast  in  the  tissue.  This  can  be  seen 
in  the  reduced  amplitude  of  the  field  at  the  geometric  focus.  Note  that  this  reduction 
does  not  stem  from  increased  absorption,  but  is  due  to  the  fact  that  the  slices  were 
taken  through  the  geometric  focus.  In  other  slices  in  the  plane  it  will  appear  that  the 
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Figure  5.5:  Acoustic  pressure  and  temperature  above  background  (37°C)  at  16/is, 
24)us,  and  32yus  for  the  ±20%  inhomogeneity  contrast  case.  Note  that  the  ±20%  and 
the  ±10%  cases  give  qualitatively  similar  results,  because  the  inhomogeneities  are  in 
the  same  locations.  Axes  labeled  in  cm. 


inhomogeneous  case  fields  are  stronger  there.  Refiection  from  the  water-tissue  inter¬ 
face  can  be  seen  at  about  x=2.75  cm.  The  reflections  from  all  three  slices  coincide, 
as  the  travel  path  between  the  source  to  the  interface  through  homogeneous  water 
is  the  same.  The  temperature  profiles  along  slices  in  x  and  y  passing  through  the 
geometric  focus  are  shown  in  Figure  5.7  at  32/us.  Again,  comparisons  of  the  overall 
heating  effect  are  difficult  to  make  based  on  slices  through  any  particular  point,  as 
the  local  effects  of  the  inhomogeneity  will  affect  the  slice  profiles. 
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Figure  5.6:  Acoustic  pressure  slices  along  x  and  y  through  the  geometric  focus  for  the 
±0%,  ±10%,  and  ±20%  inhomogeneity  contrast  cases. 


The  Overall  Effect  of  Tissue  Inhomogeneity 

It  is  difficult  to  gauge  the  effect  of  tissue  inhomogeneity  on  the  heating  by  taking 
single  slices  through  the  tissue.  Instead,  we  perform  a  numerical  integral  in  2-D  space 
over  the  tissue  half-plane  to  determine  the  average  heating  that  occurred  in  this  region 
of  interest.  The  formula  used  is 

T,,i  (5.21) 

i  3 

where  the  sum  is  taken  over  the  computational  domain.  We  find  that  Tsum  is  identical, 
independent  of  the  amount  of  inhomogeneity  contrast  present.  The  reason  for  this  is 
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Figure  5.7:  Temperature  elevation  slices  along  x  and  y  through  the  geometric  focus 
for  the  ±0%,  ±10%,  and  ±20%  inhomogeneity  contrast  cases. 


that  while  the  beam  pattern  is  distorted  in  the  region  of  interest  for  the  inhomogeneous 
runs,  the  amount  of  absorption  and  the  center  frequency  remain  the  same  as  for  the 
reference  case.  Hence,  the  amount  of  acoustic  energy  deposited  is  the  same. 

The  localized  field  distortion  in  inhomogeneous  tissue  leads  to  local  regions  of 
elevated  temperature  which  can  lead  to  higher  peak  temperatures  than  for  the  zero 
contrast  case.  Small-scale  diffraction  and  lensing  effects  are  more  pronounced  for  the 
runs  with  higher  inhomogeneity  contrast.  Figure  5.8  shows  the  peak  temperature 
found  in  the  tissue  for  each  of  the  inhomogeneity  contrast  runs  for  the  same  initial 
conditions  at  the  source.  It  is  noted  that  a  monotonic  trend  to  higher  peak  local 
temperatures  is  obtained  as  inhomogeneity  contrast  increases.  This  can  have  rami- 
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Figure  5.8:  Peak  local  temperature  rise  found  in  the  tissue  for  various  inhomogeneity 
contrast  cases,  showing  a  general  trend  to  higher  peak  local  temperatures  for  higher 
inhomogeneity  contrast. 


fications  to  dose  and  safety  studies,  where  the  maximum  pressure  and  temperature 
may  be  of  concern.  Real  soft  tissues  can  be  expected  to  exhibit  about  ±10%  inhomo¬ 
geneity  contrast  in  the  values  of  their  background  propagation  parameters  at  constant 
temperature,  depending  on  the  organ  [43]. 

5.4-3  Severe  Inhomogeneity  Effects 

Next  we  look  briefly  at  an  example  where  a  severe  mismatch  in  properties  between 
the  tissue  and  some  inclusion  in  the  tissue  exists.  This  may  occur  if  a  cavitation 
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cloud,  a  gas  pocket,  or  a  bone  cross  section  are  in  the  path  of  the  propagation.  In  our 
example  we  assume  the  properties  within  the  inclusion  are  homogeneous,  and  have 
the  values 


Figure  5.9:  Pressure  and  temperature  elevation  for  the  ca.se  where  a  severely  mis¬ 
matched  cavity  inhomogeneity  lies  in  the  tissue  in  the  vicinity  of  the  beam.  Axes 
labeled  in  cm. 

c  =  250  m/s, 
p  =  500  kg/m^ 
a  =  1  Np/m, 

(3  =  10. 


(5.22) 
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These  values  attempt  to  represent  parameters  which  might  be  found  in  a  region 
containing  cavitation  bubbles.  No  details  of  the  properties  of  such  an  inclusion  are 
available  at  the  time  of  writing,  but  it  can  be  expected  that  both  the  sound  speed 
and  the  density  of  the  effective  medium  are  lower  than  that  of  the  tissue,  and  that 
the  nonlinearity  coefficient  would  be  significantly  higher  than  that  of  the  tissue  [42]. 

Figure  5.9  shows  simulation  output  from  an  example  with  strong  scattering.  The 
severe  pc  mismatch  displaces  the  pulse  to  the  left  (24)us),  and  leaves  a  distorted  hot 
spot  (32;us).  The  tissue  inhomogeneity  contrast  for  this  run  was  ±10%.  Output  such 
as  this  suggests  that  focused  ultrasound  surgery  is  most  accurate  in  soft  tissue  with 
minimum  obstructions.  Any  significantly  mismatched  inclusion  will  move  a  typical 
size  therapeutic  beam  off  course  by  an  amount  comparable  to  the  size  of  the  inclusion. 
Further,  treatment  of  tissue  in  the  inclusion’s  shadow  is  precluded  if  the  mismatch 
leads  to  strong  scattering  of  the  sound. 

5.5  CW  Heating  in  Quiescent  Tissue 

We  now  consider  tissue  heating  due  to  continuous  wave  (CW)  insonation.  The  in- 
sonations  presented  in  the  following  sections  will  have  longer  durations  than  the  single 
pulses  seen  earlier.  For  insonations  of  severed  seconds  duration  we  can  expect  real  tis¬ 
sue  to  experience  physical  changes  in  its  material  properties.  These  time-varying 
background  (TVB)  changes  will  be  discussed  in  the  next  section,  but  for  now  we  limit 
our  study  to  quiescent,  or  unchanging  tissue.  For  these  simulations  we  use  the  same 
acoustic  and  thermal  algorithms  developed  earlier  in  this  chapter,  but  we  allow  the 
source  to  be  active  for  a  longer  duration  (100  ps)  so  that  the  entire  computational 
domain  can  be  filled  with  sound  and  reach  a  time-harmonic  steady  state.  Figure  5.10 
shows  snapshots  of  the  pressure  and  temperature  for  a  100  ps  pulse,  which  is  long 
enough  for  the  acoustic  pressure  field  to  reach  a  time-harmonic  steady  state,  as  prop¬ 
agation  across  the  computational  domain  requires  only  30  ps.  Absorbing  boundary 
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conditions  on  all  outer  edges  of  the  computational  domain  act  to  minimize  reflected 
sound  waves  into  the  interior  of  the  domain.  Since  all  the  remaining  calculations 
require  knowledge  of  the  acoustic  pressure  field  only  in  addition  to  the  tissue  prop¬ 
erties,  we  can  use  the  results  from  the  long  pulse  simulation  to  infer  the  behavior 
of  the  pressure  and  temperature  fields  for  any  time  period  desired.  This  assumption 
that  the  pressure  field  will  look  the  same  from  one  acoustic  cycle  to  the  next  is  what 
we  call  continuous  wave  (CW)  behavior  in  a  quiescent  medium.  To  confirm  that  the 


Figure  5.10:  Pressure  and  temperature  snapshots  for  a  pulse  which  is  long  enough 
so  that  the  acoustic  field  reaches  steady  state  on  the  computational  domain.  Axes 
labeled  in  cm. 


pulse  was  long  enough  to  attain  time-harmonic  steady  state  acoustics  we  plot  an  axial 
slice  through  the  pressure  field  at  50,  75,  and  100  fis  in  Figure  5.11.  The  pressure 
field  slices  are  congruent  at  these  times,  indicating  that  the  snapshots  were  taken  at 
integral  multiples  of  a  period  and  are  at  steady  state.  Figure  5.11  also  shows  tem¬ 
perature  along  the  transducer  axis  rising  with  insonation  time,  achieving  a  value  of  1 
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2.5  3 

x(cm) 


Figure  5.11:  Pressure  and  temperature  slices  along  the  transducer  axis  taken  at  50, 
75,  and  100  fis.  Pressure  is  at  steady  state,  while  temperature  continues  to  rise  with 
insonation  time. 


mK  above  background  after  100  /is.  The  linear  rate  of  temperature  rise  is  consistent 
with  the  steady  state  conditions. 

Typically  the  FDTD  time  step  size  for  acoustic  calculations  at  1  MHz  is  on  the  or¬ 
der  of  10  nanoseconds.  Thus  it  is  not  practical  to  obtain  results  for  a  10  s  insonation 
by  running  the  basic  acousto-thermal  code  until  10  s  are  reached  (this  would  re¬ 
quire  approximately  two  years  of  execution  time  using  current  computers!).  Since  the 
thermal  problem  occurs  on  a  much  slower  time  scale,  it  is  not  necessary  to  continue 
acoustic  field  calculations  once  a  time-harmonic  steady  state  is  achieved.  Coarser  time 
gridding  is  used  to  advance  the  calculations  of  the  temperature  field.  To  obtain  the 
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temperature  for  insonations  much  longer  than  the  acoustic  simulation  would  permit, 
we  integrate  the  heat  deposited  from  the  CW  acoustic  field  over  one  period,  then  we 
assume  that  this  heat  deposition  pattern  remains  unchanged  (this  assumption  will  be 
relaxed  later)  for  the  remainder  of  the  insonation.  The  temperature  field  is  calculated 
using  the  bioheat  equation  in  a  separate  FDTD  program  using  a  much  larger  time 
step  (about  10  fis)  which  can  proceed  to  calculate  results  for  insonations  of  several 
seconds  in  about  10  minutes  of  execution  time  using  current  machines.  Thus  knowl¬ 
edge  of  the  CW  pressure  field  allows  extrapolation  to  calculate  the  heating  for  much 
longer  times,  on  the  order  of  several  seconds  for  therapeutic  applications. 

The  method  for  calculating  the  CW  heat  deposition  term  consists  of  integrating 
the  heat  deposition  field  in  time  for  one  acoustic  cycle  having  period  Ta,  after  steady 
state  is  reached. 


pto+Ta 

Q(x)acou8ticcycle=  /  Q(x,t)di,  (5.23) 

Jto 

where  to  is  great  enough  to  achieve  a  time-harmonic  steady  state  acoustic  field  on  the 
computational  domain.  Subsequent  FDTD  calculations  of  the  bioheat  equation  can 
then  use  a  larger  time  step  than  that  used  in  the  pressure  calculations,  and  the  rate 
of  heat  deposition  used  in  the  bioheat  calculations  Q{x)bioheat  has  the  same  spatial 
distribution,  but  has  magnitude 


r 

^/  \  _ A  ^tbioheat 

Vv^jbioheat  —  acoustic  cycle  ? 

Ta 


(5.24) 


where  ^ibioheat  is  the  timestep  used  for  the  bioheat  equation  solver. 

The  10  s  time  trace  of  the  temperature  at  the  geometric  focus  of  the  source  is 
shown  in  Figure  5.12,  where  the  insonation  lasted  for  5  s  followed  by  a  cooling  time 
of  5  s.  The  simulation  assumes  tissue  parameters  for  soft  tissue  as  listed  in  equation 
(5.20).  The  spherical  section  bowl  transducer  emits  a  1  MHz  field  at  1  MPa  source 
pressure.  The  aperture  of  the  bowl  is  4  cm  and  the  radius  of  curvature  is  3  cm. 

At  this  point  it  is  appropriate  to  examine  the  bioheat  equation  (5.1)  again  for 
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Figure  5.12:  Temperature  at  the  geometric  focus  due  to  a  10-second  simulation  for 
a  1  MHz  bowl  source  at  1  MPa  in  soft  tissue.  The  source  was  ON  for  5  seconds, 
followed  by  5  seconds  with  the  source  OFF. 


insight  into  the  individual  terms.  We  rewrite  the  bioheat  equation  here  for  convenience 

dT  _  Ktiss  ^rp  _  ^  _j_  Qbioheat  25^ 

dt  pCtiss  P^tiss  P^tiss 

A  number  of  effects  are  represented  in  (5.25).  The  driving  source  for  the  heating 
of  the  tissue  is  Q  which  is  a  function  of  space  only  for  steady  heating  in  quiescent 
media.  The  other  terms  (from  left  to  right)  represent  temperature  rise,  conduction 
due  to  the  spatial  gradient  of  temperature,  and  perfusion  due  to  the  elevation  in 
temperature  above  the  ambient  temperature,  Ta  =  S7°C.  Figure  5.13  shows  results 
from  a  1-second  simulation  using  the  same  bowl  source  described  in  the  previous 
paragraph.  The  source  is  on  for  0.5  seconds  and  turned  off  for  the  remaining  0.5 
seconds.  The  figure  plots  the  temperature  at  the  location  of  the  transducer  focus 
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under  three  conditions: 


Figure  5.13:  The  time  trace  of  temperature  at  the  focus  for  a  0.5  s  on  then  0.5  s 
off  insonation.  The  traces  show  the  effect  of  conduction  jSTtiss  and  perfusion  Wb  on 
temperature. 


1.  Kt\ss  =  Wb  =  0:  In  this  case  no  conduction  or  perfusion  axe  simulated,  and 
the  only  effect  is  due  to  heating  by  the  source  Q.  The  bioheat  equation  is  thus 
reduced  to  the  linear  relation 


dT 

dt 


Q 


pC^. 


tiss 


(5.26) 


We  see  for  this  case  that  the  temperature  rise  is  linear,  with  a  slope  of  QI{pCt\s^, 
and  that  no  cooling  occurs  after  the  source  is  extinguished.  In  other  words,  if 
integrated  once  over  time  r 

Q  . 


T(r)  =  7;  + 


pQ 


tiss 


(5.27) 
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Figure  5.14:  Temperature  difference  between  the  focal  temperature  with  and  without 
considering  perfusion.  An  inflection  point  in  the  curve  at  0.5  s  indicates  where  (T—Ta) 
stops  increasing  and  starts  decreasing. 


2.  Ktiss  >  0,  Wfc  =  0:  If  conduction  is  taken  into  account,  but  not  perfusion,  we 
see  a  nonlinear  profile  to  the  temperature  of  the  hot  spot  in  Figure  5.13.  The 
conduction  term  acts  to  remove  some  of  the  heat  from  the  immediate  vicinity 
of  the  focus  during  and  after  the  insonation.  ATtiss  used  was  0.6W/(m  •  K). 

3.  Ktiss  >  0,1^6  >  0:  For  the  case  where  both  conduction  and  perfusion  are 

taken  into  account,  we  see  a  curve  almost  identical  to  that  for  the  iTtiss  > 
0,  =  0  case.  However,  to  look  at  whether  the  perfusion  has  any  effect  on 
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the  temperature  profile  we  plotted  the  difference  between  the  focal  temperature 
with  and  without  perfusion.  The  result  is  shown  in  Figure  5.14  and  shows  that 
indeed  there  is  a  small  effect  attributable  to  perfusion.  Figure  5.14  indicates 
that  the  perfusion  effect  becomes  increasingly  important  as  the  temperature 
rises  above  ambient  temperature,  because  the  perfusion  term  contains  the  factor 
{T  —  Ta).  Wb  results  from  heat  exchange  from  the  heated  tissue  to  the  ambient 
blood  flow  network,  assumed  uniform  in  the  tissue.  The  value  used  for  Wb  was 
0.5kgl{m^  ■  s). 

5.6  Therapeutic  Ultrasound  in  Nonlinear  Time-Varying  Tissue:  CW 
Simulations 

A  model  for  propagation  of  finite  amplitude  ultrasonic  waves  in  nonlinear,  absorbing, 
time-varying  background  (TVB)  media  was  presented  in  Chapter  2.  Numerical  sim¬ 
ulations  using  the  FDTD  method  are  used  to  investigate  the  thermal  and  acoustic 
effects  from  focused  transducers  in  media  with  TVB  sound  speed  and  attenuation 
coefficient. 

The  motivation  for  these  simulations  is  the  experimental  data  showing  the  temperature- 
dependence  of  Co  and  ao  [4,  17].  The  data  measured  by  Bamber  and  Hill  and  by 
Damianou,  et  al.  are  adapted  for  the  present  study  by  using  one  data  point  per  5°C 
between  30°  C  and  90°  C.  The  results  are  shown  in  Figure  5.15  along  with  polynomial 
fits  to  the  data.  The  polynomial  used  to  fit  the  sound  speed  to  the  data  is 

c{T)  =  1466  4-  8.43r  -  O.mT^  -t- 1.27  x  lO^^T^  -  3.33  x  lO^^T^  m/s.  (5.28) 

The  polynomial  used  to  fit  the  absorption  coefficient  to  the  data  is 

a(T)  =  -43.8  +  6.24T  -  0.303r  -|-  7.06  x  10“^T^  -  8.39  x  10~®r^ 

^  ^  (5.29) 

-b  4.91  X  10“^r^  -  1.12  X  10-®T®  Np/m. 
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Figure  5.15;  Temperature  dependence  of  sound  speed  and  absorption  coefficient  in 
soft  tissue  taken  from  published  laboratory  measurements  (symbols)  and  the  corre¬ 
sponding  polynomial  fits  to  the  data  (solid). 


Caution  should  be  exercised  not  to  extrapolate  the  polynomials  outside  the  region 
of  fitting,  as  they  only  model  the  data  inside  the  interpolation  regions,  here  30° C  to 
90°  C. 

These  data  imply  that  during  the  course  of  a  single  focused  ultrasound  treatment 
tissue  properties  will  change  and  quite  significantly  for  the  absorption  coefficient  Oq  if 
temperature  exceeds  60° C.  Data  for  other  properties  such  as  nonlinearity  coefficient 
and  density  are  not  available.  The  temperature-dependence  of  oq  implies  the  rate  of 
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heat  deposition  will  also  be  temperature  (and  time)  dependent,  since  the  heating  rate 
is  directly  proportional  to  the  absorption  coefficient. 

To  estimate  the  importance  of  the  TVB  behavior  on  the  acoustic  waves,  we  used 
order  of  magnitude  estimates  in  section  2.5.2  to  obtain  the  relative  magnitudes  of 
the  TVB  terms  in  the  TVB  wave  equation  (2.56)  compared  to  the  other  terms  in 
the  wave  equation.  For  tissue  insonated  by  a  1  MHz  bowl  source  of  2  cm  aperture, 
having  a  source  pressure  of  1  MPa,  we  concluded  that  for  the  range  of  parameters 
of  interest,  the  time  variation  in  the  tissue  background  properties  is  not  important 
at  the  acoustic  time  scales.  In  other  words,  we  formally  demonstrated  that  the  slow 
TVB  effects  (temperature  dependence  with  time  scales  of  1  s)  are  too  slow  to  affect 
the  acoustic  waves  having  time  scales  on  the  order  of  1  fis. 

One  may  be  inclined  to  conclude  that  the  previous  paragraph  implies  that  the 
TVB  is  not  a  factor  in  acoustic  applications  of  ultrasound.  This  is  however  untrue, 
since  a  typical  treatment  does  occur  over  a  span  of  several  seconds  at  the  least,  to 
achieve  a  useful  temperature  rise.  Thus  the  variation  in  the  background  acoustic 
properties  of  the  tissue  should  be  taken  into  consideration  for  insonations  with  dura¬ 
tions  comparable  to  the  thermal  time  scales.  In  this  chapter  we  will  show  the  results 
of  simulations  which  couple  the  acoustic  and  the  thermal  behavior  of  tissue  to  see 
what  the  effect  on  overall  heating  is. 

5.6.1  An  Example  of  TVB  Behavior 

Now  we  look  at  one  case  where  the  background  sound  speed  and  attenuation  coef¬ 
ficient  are  allowed  to  vary  in  time  according  to  the  temperature  of  the  tissue  as  a 
function  of  space  and  time.  In  this  case  the  sound  field  can  be  expected  to  change 
with  both  the  thermal  time  scale  as  well  as  the  acoustic  time  scale.  Specifically,  as 
temperature  increases  in  the  focal  region,  we  can  expect  the  rate  of  temperature  rise 
to  increase  because  the  driving  term  in  the  bioheat  equation  is  dependent  on  a,  which 
is  monotonically  increasing  with  temperature.  In  fact,  a  doubling  of  the  attenuation 
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coefficient  appears  to  be  possible  over  the  period  of  a  hyperthermia  treatment  of  sev¬ 
eral  seconds  duration.  This  is  because  the  temperature  range  up  to  70°  C  and  higher 
is  possible  with  an  insonation  of  a  few  seconds  from  a  typical  laboratory  hyperthermia 
focused  source. 


Figure  5.16:  Flowchart  showing  the  iterative  method  for  coupling  the  pressure  and 
temperature  calculations  in  the  TVB  simulations. 


The  flowchart  in  Figure  5.16  shows  how  the  acoustic  and  the  thermal  solvers  were 
coupled  for  the  TVB  computations  via  the  CW  heating  term  Q  as  before.  Periodically 
(at  1  s  intervals)  the  thermal  solver  is  made  to  output  updated  background  properties 
of  the  tissue  to  data  files.  The  sound  speed  and  attenuation  coefficient  are  updated 
using  the  polynomial  fitting  routine  derived  from  the  laboratory  measurements  of  the 
temperature-dependence  of  c  and  a.  The  updated  c  and  a  are  then  used  as  inputs 
to  the  acoustic  solver  to  calculate  a  new  pressure  field,  and  a  new  Q  field,  and  so  on 
in  an  iterative  fashion.  Thus  the  pressure  field,  which  drives  the  heating,  is  in  turn 
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Figure  5.17:  Temperature  profiles  across  the  axis  in  the  focal  plane  at  3,  4,  5  s, 
showing  increased  heating  in  the  TVB  case. 


affected  by  the  temperature  field.  The  temperature’s  effect  on  a  can  be  particularly 
important  in  hyperthermia  applications. 

In  Figure  5.17  we  show  slices  across  the  axis  of  a  bowl  transducer  as  they  evolve 
in  time  for  quiescent  and  TVB  simulations.  Note  that  at  3  s  the  TVB  case  results 
are  close  to  the  quiescent  results,  or  even  a  little  lower  in  temperature.  The  reason 
is  that  at  3  s  the  temperature  is  about  55°  C,  and  at  this  temperature  the  absorption 
coefficient  has  dropped  slightly  compared  to  that  at  37°  C.  As  the  insonation  time 
increases,  the  temperature  exceeds  60°C,  where  there  is  a  dramatic  increase  in  a. 
In  this  case  the  temperature  from  the  TVB  simulation  is  higher  than  that  from  the 
quiescent  simulation.  This  result  can  be  explained  if  we  look  at  the  profile  of  the 
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Figure  5.18:  Slices  of  Q  taken  across  the  axis  at  the  focal  zone  at  1  s  and  5  s  for  the 
TVB  case.  A  dramatic  increase  in  the  heat  deposition  rate  Q  can  be  observed  if  the 
attenuation  coefficient  (directly  proportional  to  Q)  is  included  in  the  simulation.  Q 
would  be  steady  for  quiescent  simulations. 


driving  source  term  Q  in  the  bioheat  equation  as  it  evolves  in  the  TVB  case.  Figure 
5.18  shows  the  Q  profile  across  the  transducer  axis  in  the  focal  plane  at  1  s  and  5 
s.  Q  increases  due  to  the  increase  in  oc  with  temperature.  The  direct  proportionality 
between  a  and  Q  is  responsible  for  the  added  heating  noted  at  the  focus. 


5.7  Conclusions 

This  chapter  explored  several  facets  of  focused  ultrasound  surgery  in  soft  tissues  using 
modeling  and  numerical  simulations.  The  techniques  developed  in  earlier  chapters  for 
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modeling  acoustic  sources  and  fields  was  coupled  to  the  bioheat  equation  through  the 
acoustic  energy  deposition  driving  term  in  the  bioheat  equation.  Various  effects  were 
studied  including  the  effect  of  finite-amplitude  propagation,  inhomogeneity  contrast, 
and  slow  TVB  variation  in  temperature-dependent  tissue. 

We  showed  that  finite-amplitude  propagation  leads  to  enhanced  heating.  Be¬ 
cause  absorption  in  lossy  media  is  monotonic  with  frequency,  a  waveform  containing 
higher-frequency  harmonics  is  absorbed  more  readily  than  one  containing  only  the 
lower-frequency  fundamental.  This  result  is  not  new,  but  was  confirmed  here  for 
focused  sources  in  tissue-like  media  using  the  FDTD  simulations.  The  amount  of 
excess  heating  in  finite-amplitude  situations  over  linear  ones  depends  on  the  same 
factors  governing  the  generation  of  harmonics:  propagation  distance  to  target,  source 
frequency  and  amplitude,  and  nonlinearity  coefficient  of  the  tissue.  For  the  example 
used  in  this  study,  the  nonlinear  temperature  elevation  was  about  up  to  1.8  times 
that  of  the  linear  case  for  similar  insonations  from  a  short  pulse. 

Inhomogeneity  contrast  about  a  mean  value  wets  studied  because  of  its  effect  on 
the  profile  of  the  acoustic  beam  and  the  resulting  distortion  of  the  hot  spot  near  the 
focus.  The  present  study  used  a  tissue  inhomogeneity  model  derived  from  laboratory 
measurements  using  human  tissue  samples  [43].  It  was  found  that  while  increasing  the 
percent  inhomogeneity  contrast  does  not  affect  the  overall  thermal  energy  deposition 
in  the  extended  tissue  sample  (a  statement  of  conservation  of  energy).  However,  the 
peak  temperatures  in  the  computational  domain  increased  monotonically  with  percent 
inhomogeneity  contrast.  We  hypothesize  that  the  cause  for  this  is  local  lensing  of  the 
acoustic  field,  and  accompanying  local  heating  due  to  the  inhomogeneity  result  when 
small-scale  inhomogeneities  are  inserted  in  the  propagation  path.  The  peak  local 
temperature  excursions  can  have  implications  for  safety  studies  of  FUS  devices.  Severe 
inhomogeneities,  such  as  obstruction  by  bone,  vasculature,  bladders,  or  bubble  clouds 
can  have  dramatic  effects  on  the  acoustic  and  temperature  fields  due  to  reflections 
and  scattering  from  the  obstruction. 
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Finally,  when  the  effect  of  temperature-dependence  on  tissue  sound  speed  and 
attenuation  coefficients  was  taken  into  account  (TVB  simulations),  increased  heat¬ 
ing  near  the  focus  was  observed.  It  has  been  observed  in  the  laboratory  that  the 
attenuation  coefficient  of  soft  tissue  is  approximately  monotonically  dependent  on 
temperature  [17]  between  SO^C  and  90°C,  with  the  absorption  coefficient  nearly  dou¬ 
bling  between  50°  C  and  70°  C.  This  effect  will  cause  an  increase  in  the  rate  of  heating 
in  the  regions  near  the  focus  of  a  FUS  device.  The  present  study  shows  that  since  the 
heating  rate  is  directly  proportional  to  the  absorption  coefficient,  slow  TVB  effects 
need  to  be  accounted  for  in  modeling  and  simulating  FUS  treatments. 
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Chapter  6 

CONCLUSIONS 

This  dissertation  addressed  issues  associated  with  forming  a  high-amplitude  focus 
using  arrays  and  extended  focused  sources  in  various  propagation  media.  Specifically, 
we  addressed  the  problems  which  may  face  time  reversal  focusing  in  the  presence  of 
debilitating  effects  such  as  phase  jitter,  absorption,  and  nonlinearity.  Further,  the 
behavior  of  focused  sources  used  for  tissue  heating  in  biomedical  applications  was 
examined  in  the  presence  of  nonlinearity,  inhomogeneity,  and  time-varying  medium 
properties.  These  studies  were  conducted  under  the  umbrella  of  finite-amplitude 
acoustic  models  and  associated  numerical  simulation  codes. 

The  practical  applications  of  acoustics  and  ultrasonics  in  medical,  industrial  and 
military  uses  are  almost  completely  restricted  to  two  classes  of  applications: 

1.  Nondestructive  diagnostic  and  imaging  applications.  Diagnostic  medical  imag¬ 
ing  with  ultrasound,  nondestructive  testing  of  materials  (NDT),  and  communi¬ 
cations  and  surveillance  are  examples  of  this  use. 

2.  Applications  requiring  focused,  high-intensity  sound  for  the  purpose  of  affecting 
physical  change  to  the  objects  insonated  near  the  focal  spot.  Examples  are  fo¬ 
cused  ultrasonic  surgery,  ultrasonic  welding  of  plastics,  and  mine  neutralization. 

This  dissertation  deals  with  the  latter  use  for  acoustic  arrays  and  focused  sources 
in  underwater  and  biomedical  applications.  Several  topics  in  linear  and  nonlinear 
acoustic  focusing  were  studied,  with  emphasis  on  the  evaluation  of  the  focusing  char¬ 
acteristics  of  time  reversal  arrays  (TRA’s)  and  therapeutic  hyperthermia  sources  for 
focused  ultrasound  surgery  (FUS). 
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Two  recently-proposed  applications  for  acoustic  arrays  and  focused  sources  are 
high-intensity  time  reversal  arrays  and  focused  hyperthermia  applicators.  The  un¬ 
derlying  concepts  for  the  operation  of  these  devices  are  reasonably  well  understood. 
However,  both  of  these  technologies  in  their  current  forms  are  relatively  recent  addi¬ 
tions  to  the  arsenal  of  underwater  and  biomedical  engineers,  and  until  now  they  have 
not  become  commercially  viable.  We  now  examine  the  background  and  the  contri¬ 
butions  of  this  dissertation  to  each  of  the  main  problems  addressed  in  this  study,  as 
well  as  future  directions  for  the  use  of  TRA’s  and  FUS  sources. 

6.1  Time  Reversal  Arrays 

6.1.1  Motivation  for  the  TR A  Study 

The  advent  of  digital  data  aquisition  and  computer-controlled  transducer  arrays  al¬ 
lowed  scientists  to  realize  the  acoustic  phase  conjugate  mirror,  or  time  reversal  aj-ray. 
The  device  allows  for  automatic  phase  correction  at  the  array  to  produce  a  focus 
at  the  location  of  a  scatterer  in  an  unknown  inhomogeneous  medium  containing  ar¬ 
bitrary  multipath  and  multiple  scattering.  The  technique  relies  on  the  stationary 
medium  providing  a  matched  filter  for  the  propagation  of  acoustic  waves  from  and 
to  an  illuminated  scatterer,  even  if  the  nature  of  the  propagation  path  is  not  known. 
The  technique  is  superior  to  conventional  time-of-flight  initial  delay  focusing  because 
it  captures  the  initial  delay  information  as  well  as  detailed  waveform  distortion  in¬ 
formation  from  all  propagation  paths,  and  corrects  for  both  using  the  medium  as  a 
spatial-temporal  filter.  The  technique  was  hailed  as  a  true  advance  in  imaging  and 
focusing  applications  because  it  permitted  the  phase  correction  due  to  unknown  ex¬ 
tended  inhomogeneities  rather  than  just  a  thin  phase  aberrating  layer  near  the  array. 
One  intriguing  aspect  of  time  reversal  array  technology  is  the  promise  of  devices  that 
can  direct  intense  acoustic  energy  onto  one  or  more  scattering  targets.  The  idea  is  to 
deposit  sufficiently  intense  fields  onto  the  targets  to  cause  permanent  physical  change 
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or  destruction.  This  idea  could  be  useful  for  accurate  targeting  and  destruction  of 
kidney  stones,  or  for  remote  mine  neutralization  at  sea. 

Currently,  the  most  successful  means  for  removal  of  kidney  stones  involves  imaging 
with  X-rays,  then  placing  the  patient  and  the  stone  in  a  coupling  water  bath  at  the 
focus  of  a  lithotripter  source.  The  lithotripter  machine  then  fires  hundreds  of  shots 
(focused  shock  waves)  at  the  stone.  While  lithotripsy  is  effective  and  is  in  wide  use, 
several  problems  plague  lithotripters  and  reduce  their  effectiveness.  The  targeting 
requires  the  use  of  X-ray  radiation  to  locate  the  stones.  In  addition,  the  targeting 
assumes  that  the  propagation  path  is  a  homogeneous  medium,  and  thus  both  the 
position  of  the  focus  and  the  shape  of  the  beam  will  suffer  from  variations  in  the 
index  of  refraction  resulting  from  the  propagation  through  the  human  body  and  the 
body-water  interface.  Finally,  the  stone  is  assumed  to  remain  stationary  for  the 
duration  of  the  treatment,  or  for  a  portion  of  the  treatment  between  repositioning.  It 
is  thought  that  a  large  percentage  of  shots  fired  at  a  kidney  stone  miss  their  intended 
target  due  to  breathing  motion  and  other  targeting  errors  [67]. 

Using  time  reversal  arrays  has  been  suggested  as  a  means  of  targeting  kidney 
stones  that  can  correct  for  the  unknown  refractive  index  variations  in  the  patient’s 
body.  Iterative  use  of  the  time  reversal  method  can  track  moving  scatterers  and 
maintain  the  focus  of  the  array  on  a  kidney  stone  in  a  breathing  patient.  The  next 
step  in  the  process  involves  amplifying  the  array  outputs  to  send  high-amplitude 
pulses,  or  shock  waves,  toward  the  kidney  stone  for  the  destruction  phase.  This  is 
repeated,  always  targeting  the  brightest  scatterer  until  the  procedure  is  complete. 
The  concept  has  been  successfully  demonstrated  in  the  laboratory  for  low  amplitudes 
[67].  A  similar  technique  was  proposed  for  the  remote  neutralization  of  mines  in  a 
shallow  water  channel,  and  the  proposed  solution  was  similar  to  that  intended  to 
destroy  kidney  stones. 

Ideally,  time  reversal  systems  will  operate  very  effectively  to  focus  onto  their  tar¬ 
gets  according  to  theory.  However,  in  real  applications  we  will  encounter  problems 
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that  will  degrade  the  effectiveness  of  the  TRA’s  as  focusing  implements.  These  factors 
that  compromise  a  time  reversal  system’s  performance  are  called  debilitating  factors 
in  this  dissertation.  Several  debilitating  factors  are  considered  in  this  study.  These 
include  imperfect  initial  phasing  of  the  TRA,  absorption,  and  nonlinear  propagation. 

The  arrival  times  for  different  propagation  paths  in  a  time  reversal  channel  can 
be  adversely  affected  if  jitter  exists  in  the  recorded  or  the  retransmitted  signals  from 
the  array  elements.  This  study  looks  at  the  quantitative  effect  of  this  jitter  as  a 
debilitating  factor. 

Absorption  in  the  propagation  medium  will  degrade  the  performance  of  time  re¬ 
versal  focusing  systems.  This  is  because  the  absorption  term  in  the  wave  equation 
contains  odd-order  time  derivatives,  which  are  not  invariant  under  time  reversal.  This 
has  unfortunate  implications  for  high-intensity  time  reversal  systems,  and  thus  absorp¬ 
tion  is  studied  as  a  debilitating  factor.  This  is  especially  important  for  high-amplitude 
applications. 

In  order  for  high-amplitude  time  reversal  systems  to  be  realized,  we  must  re¬ 
examine  the  assumptions  and  the  mathematics  of  the  phase-conjugate  arrays  with 
finite-amplitude  acoustics  taken  into  consideration.  As  source  amplitude  and  propa¬ 
gation  distance  increase,  energy  is  transfered  from  the  lower  frequency  fundamental 
into  higher  frequency  harmonics  due  to  nonlinearity.  In  fluids  the  energy  loss  is  invari¬ 
ably  higher  for  higher  frequencies.  The  result  is  an  overall  loss  of  focal  spot  pressure 
amplitude  when  finite-amplitude  time  reversal  is  used  in  lossy  media.  For  applica¬ 
tions  requiring  high  acoustic  pressures  at  the  focus  like  those  described  above,  this 
means  that  a  finite-amplitude  time  reversal  system  will  not  be  as  efficient  as  a  lin¬ 
ear  one.  This  dissertation  explores  this  degradation  in  a  quantitative  way,  and  uses 
simulations  to  provide  evidence  of  the  extent  of  the  degradation  as  a  result  of  the 
debilitating  factor:  finite-amplitude  propagation. 

In  order  to  get  an  idea  of  how  well  we  may  expect  real  time  reversal  systems  to 
perform  for  underwater  and  biomedical  applications,  this  dissertation  examined  the 
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operation  of  time  reversal  arrays  under  nonideal  conditions,  i.e.  in  the  presence  of 
the  debilitating  effects  above. 

6.1.2  Contributions  Made  to  the  TRA  Problem 

Errors  in  the  initial  phasing  of  time  reversed  signals  can  occur  if  the  electronics  of  the 
time  reversal  system  (aquisition,  storage,  signal  processing)  are  of  poor  resolution. 
Depending  on  the  frequency  and  bandwidth  of  the  signals,  errors  in  phasing  at  the 
array  can  result  in  loss  of  focusing  by  a  TRA.  Jitter  in  the  retransmitted  phases  of  the 
time  reversed  signals  was  found  to  be  tolerable  up  to  about  one-eighth  to  one-sixth  of 
a  fundamental  period,  beyond  which  a  loss  of  focus  occurs  [41].  This  is  corroborated 
by  studies  of  time-delay  phased  arrays  by  other  researchers  [70,  62],  but  the  present 
study  looked  explicitly  at  the  effect  of  jitter  in  time  reversal  applications. 

Absorption  in  the  propagation  medium  is  a  serious  factor  in  degrading  the  ability 
of  TRA’s  to  place  an  intense  focus  at  the  desired  target  location.  This  is  only  par¬ 
tially  due  to  the  loss  of  acoustic  energy  in  the  linear  absorption  in  the  thermoviscous 
medium.  What  makes  absorption  a  debilitating  factor  for  TRA’s  is  the  fact  that  the 
absorption  term  in  the  wave  equation  causes  a  violation  of  time  reversal  invariance. 
The  wave  equation  will  no  longer  have  twin  solutions  (forward  time  and  backwaxd 
time  solutions).  As  a  result,  the  overall  ability  of  a  TRA  to  form  a  sharp  intense  focus 
is  reduced.  This  has  been  known  to  be  the  case  in  absor]}ing  media,  but  this  study 
presented  formal  analysis  of  the  wave  equation  for  thermoviscous  fluids  showing  the 
cause  for  the  time  reversal  invariance  violation,  along  with  examples  of  this  effect  in 
simulations  of  TRA’s  for  biomedical  applications. 

While  absorption  is  found  to  be  a  debilitating  factor  for  time  reversal  arrays.  This 
study  used  propagation  through  a  model  section  of  the  human  skull  and  brain  as  a 
tool  to  examine  the  effects  of  absorption  in  the  skull  bone  layer  on  a  TRA  placed  in 
water  outside  the  subject’s  head.  It  was  found  that  for  the  pulse  used  in  this  study  a 
marked  reduction  in  focal  amplitude  and  an  increase  in  full-width-at-half-maximum 
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occurred. 

The  operation  of  a  time  reversal  array  was  re-examined  using  the  nonlinear  ab¬ 
sorbing  wave  equation  as  a  model  for  wave  propagation.  It  was  shown  that  the  non¬ 
linearity  in  itself  is  not  responsible  for  degrading  the  performance  of  TRA  systems. 
Only  when  combined  with  absorption  in  the  medium  or  amplification  at  the  array 
were  finite-amplitude  TRA  systems  adversely  affected.  It  was  shown  analytically 
that  the  extra  absorption  loss  due  to  high  frequency  content  of  shocked  waveforms 
in  nonlinear  TRA’s  can  have  a  very  detrimental  effect  on  the  efficiency  of  a  focusing 
system  using  TRA’s.  A  source-target  separation  of  more  than  one  to  two  shock  for¬ 
mation  distances  was  found  to  lead  to  a  reduction  of  the  peak  pressure  amplitude  at 
a  TRA  focus  to  a  fraction  of  that  expected  from  linear  propagation  alone  in  a  similar 
linear  medium.  Furthermore,  if  the  amplitude  of  the  received  signal  was  altered  at  the 
array  (amplified  for  example)  nonlinear  TRA’s  would  suffer  a  performance  penalty 
because  the  exact  replication  of  the  received  waveform  (time-reversed)  is  a  condition 
for  perfect  phase  conjugation.  This  amplification  at  the  array,  which  is  commonly 
assumed  for  proposed  applications  in  underwater  and  biomedical  applications,  is  a 
violation  of  time  reversal  invariance  and  is  a  debilitating  effect. 

6.2  Therapeutic  Ultrasound:  Focused  Ultrasound  Surgery 

6.2.1  Motivation  for  the  FUS  Study 

The  use  of  focused  ultrasound  sources  to  deposit  thermal  energy  for  therapeutic 
biomedical  applications  (hyperthermia)  was  also  investigated  in  this  dissertation. 
Much  promise  for  treatment  of  deep-seated  tumors  has  been  evidenced  by  in  vitro 
and  in  vivo  experiments  in  the  last  decade  using  focused  transducers  operating  in  the 
megahertz  frequency  range  [47,  65,  63].  The  absorbing  tissue  acts  to  convert  acoustic 
energy  to  thermal  energy,  inducing  a  temperature  rise.  The  temperature  rise  near 
the  focus  can  kill  living  cells  via  necrosis  by  exceeding  the  thermal  dose  threshold 
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associated  with  cell  death.  While  the  basic  operation  of  hyperthermia  devices  is  un¬ 
derstood,  and  models  exist  to  predict  the  behavior  of  the  acoustic  and  thermal  fields 
of  such  devices,  some  aspects  of  high-intensity  focused  surgery  systems  still  require 
further  study  in  order  to  be  made  into  viable  medical  systems.  This  dissertation 
examines  the  behavior  of  focused  ultrasonic  surgery  systems  through  numerical  simu¬ 
lations  and  modeling.  Key  aspects  of  ultrasound  focused  surgery  such  as  the  heating 
from  finite-amplitude  sources,  the  heating  in  inhomogeneous  tissue-like  media,  and 
the  heating  of  time-varying  tissue  were  studied. 

It  is  known  that  increased  tissue  heating  will  occur  from  high-amplitude  focused 
sources  due  to  the  higher  harmonic  content  of  the  acoustic  field  near  the  focus.  The 
higher  frequency  content  in  shocked  waveforms  deposit  more  energy  in  the  absorbing 
tissue  because  of  the  frequency  dependence  of  the  absorption  term  in  the  wave  equa¬ 
tion.  Predictions  of  heating  based  only  on  linear  field  calculations  in  homogeneous 
media  are  common  in  hyperthermia  research,  and  a  good  balance  between  clinical  util¬ 
ity  and  well-founded  theoretical  physical  acoustics  is  required  for  the  next-generation 
simulations  of  hyperthermia  systems. 

One  aspect  of  propagation  in  real  tissues  that  plays  an  important  role  in  hyper¬ 
thermia  is  the  presence  of  inhomogeneities  in  the  background  properties  of  soft  tissue 
and  other  organs  in  the  body.  Spatial  variation  in  the  propagation  medium  properties 
alter  the  exact  shape  and  location  of  the  acoustic  and  thermal  hot  spots.  Most  cur¬ 
rent  hyperthermia  research  is  conducted  under  controlled  conditions  using  soft  tissue 
where  the  effects  of  propagation  uncertainties  are  minimized.  Realtime  measurement 
of  such  inhomogeneities  for  case-specific  correction  would  require  sophisticated  imag¬ 
ing  techniques  in  order  to  faithfully  capture  the  nature  of  the  inhomogeneity,  as  well 
as  powerful  computational  abilities.  This  is  especially  important  if  passing  acoustic 
beams  through  regions  of  high  inhomogeneity  contrast,  or  severe  inhomogeneities  such 
as  bone,  vasculature,  bladders,  or  gas  pockets,  in  which  case  the  passage  of  sound  may 
be  obstructed  outright  and  scattered  to  unexpected  locations. 
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Another  important  feature  of  tissues  is  the  way  that  they  respond  to  heating. 
Perfusion,  the  adaptive  cooling  of  bulk  tissue  by  blood  flow,  is  temperature  depen¬ 
dent,  and  makes  for  a  challenging  control  problem.  More  important  is  the  behavior 
of  temperature-dependent  background  propagation  properties.  The  sound  speed  and 
absorption  coefficient  have  been  measured  to  be  a  function  of  temperature,  so  the 
position  and  shape  of  the  focus  will  vary  in  time  as  the  tissue  is  heated.  For  accurate 
simulations  of  bioacoustics  problems,  other  acoustic  and  thermal  properties’  charac¬ 
teristics  should  be  measured  to  be  included  in  the  calculations.  At  this  time  most 
properties  have  only  been  reported  for  body  temperature  conditions. 

6.2.2  Contributions  Mode  to  the  FUS  Problem 

Simulations  of  the  spatial  and  temporal  evolution  of  the  pressure  and  temperature 
fields  from  hyperthermia  sources  were  presented.  The  calculated  pressure  field  was 
then  used  to  determine  the  local  heating  due  to  absorption.  This  heating  was  used 
as  a  source  term  in  solutions  of  the  transient  bioheat  equation.  Simulations  for  finite- 
amplitude  FUS  sources  were  performed,  and  the  results  compared  to  similar  sources 
in  a  linear  medium.  Enhanced  heating  depended  on  the  source  pressure,  amplitude, 
frequency,  and  geometry,  as  the  generation  of  harmonics  depends  on  these  factors.  In 
simulations  performed  for  this  research  the  heating  could  be  increased  by  a  factor  of 
1.5  to  3  in  comparison  to  predictions  by  linear  models. 

We  demonstrated  that  inhomogeneities  in  the  propagation  medium’s  background 
parameters  lead  to  local  distortion  of  the  acoustic  beam,  and  a  resulting  distortion 
of  the  heating  pattern  calculated  for  the  homogeneous  (zero  contrast)  case.  The  sim¬ 
ulations  used  spatial  distribution  data  for  real  human  tissue  properties  found  in  the 
literature.  This  study  found  that  the  overall  heat  deposited  by  a  beam  in  an  inho¬ 
mogeneous  medium  with  zero-mean-plus-average  inhomogeneities  is  the  same  for  any 
amount  of  inhomogeneity  contrast.  This  is  a  statement  of  conservation  of  energy,  since 
the  acoustic  energy  loss  to  heat  will  be  the  same,  but  merely  redistributed  in  space. 
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The  peak  temperatures  observed  however  monotonically  increased  with  inhomogene¬ 
ity  contrast.  That  is  to  say,  as  inhomogeneity  contrast  was  increased,  higher  and 
higher  values  of  peak  local  temperature  were  observed  in  the  computational  domain. 
Given  the  previously-mentioned  conservation  statement,  this  implies  that  other  loca¬ 
tions  underwent  less  heating.  This  finding  has  implications  for  safety  studies,  where 
the  peak  temperature  excursions  are  of  concern,  as  well  as  for  applications  where  rais¬ 
ing  the  temperature  sufficiently  to  ensure  complete  necrosis  is  required  at  all  locations 
in  the  region  of  interest. 

Finally,  the  coupling  of  pressure  and  temperature  calculations  was  carried  out  with 
a  slow-scale  TVB  tissue  simulation.  The  tissue  was  assumed  to  change  its  background 
sound  speed  and  attenuation  coefficients  in  accordance  with  published  measurements 
for  these  quantities.  A  polynomial  fit  to  the  reported  data  was  used  to  periodically 
update  the  initial  condition  data  files  containing  the  sound  speed  and  the  attenua¬ 
tion.  This  new  data  was  then  fed  into  the  pressure  solver  program  to  compute  an 
updated  acoustic  field.  This  interdependence  of  heating  on  acoustic  pressure  and 
acoustic  pressure  on  temperature  formed  a  feedback  loop  that  dramatically  altered 
the  temperature  predictions  for  a  hyperthermia  system.  The  fact  that  the  attenuation 
coefficient  had  a  monotonically  increasing  dependence  on  temperature,  and  the  fact 
that  heating  rate  is  directly  proportional  to  absorption  coefficient  implies  an  accel¬ 
eration  in  the  rate  of  temperature  rise  for  a  system  with  slow  TVB  feedback.  The 
absorption  coefficient  can  double  over  the  span  of  a  hyperthermia  treatment  session, 
thus  doubling  the  heating  rate  for  the  same  acoustic  field.  Results  confirmed  that  as 
temperature  increased  near  the  hot  spot,  the  disparity  between  the  cases  where  TVB 
was  and  was  not  taken  into  account  grew  in  magnitude. 


116 


6.3  Other  Contributions 

This  dissertation  has  made  original  contributions  to  the  study  of  finite-amplitude 
propagation  in  inhomogeneous  thermoviscous  time-varying  background  (TVB)  media. 
Starting  from  the  equations  of  fluid  mechanics,  the  nonlinear  absorbing  wave  equation 
in  inhomogeneous  thermoviscous  fluids  was  derived  for  the  case  having  time- varying 
background  (TVB)  parameters.  This  wave  equation  was  studied  in  depth  and  nondi- 
mensionalized  to  elucidate  the  significance  of  the  TVB,  and  the  regimes  of  operation 
which  would  require  consideration  of  the  TVB  wave  equation. 

Simulations  were  carried  out  using  original  computer  codes  capable  of  representing 
the  physics  appearing  in  the  wave  equation.  Finite-difference  time-domain  (FDTD) 
simulations  allowed  detailed  study  of  the  spatial  and  temporal  behavior  of  the  sys¬ 
tems  under  consideration,  and  were  especially  helpful  in  capturing  the  essence  of  the 
behavior  of  the  focal  zone. 

6.4  Summary 

Despite  having  a  good  foundation  for  understanding  the  workings  of  time  reversal 
systems  and  acoustic  hyperthermia  processes,  this  research  reports  on  some  important 
new  findings.  In  the  case  of  time  reversal  arrays,  nonlinear  effects  were  shown  to  have 
a  serious  debilitating  effect  on  the  ability  of  a  TRA  to  form  an  intense  focus  at  the 
desired  location.  Absorption  and  amplification  of  the  time  reversed  signal  at  the 
array  were  shown  through  analysis  and  simulation  to  be  detrimental  to  the  TRA 
performance.  Jitter  in  the  time  reversed  signal  was  also  shown  to  gradually  erode  the 
sharpness  of  the  focus  of  a  TRA  system  in  an  underwater  channel.  Future  research 
involving  laboratory  measurements  for  finite-amplitude  time  reversal  systems  would 
complement  these  findings,  and  given  the  limitations  of  the  numerical  tools  available, 
perhaps  the  measurements  could  extend  the  results  to  higher  Gol’berg  numbers  and 
for  greater  distances  than  those  simulated  in  this  study.  The  simplest  form  of  a  TRA, 
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a  one-dimensional  two-element  case,  was  used  for  this  study.  It  would  be  informative 
to  use  real  multi-dimensional  arrays  in  a  laboratory  study.  Time  reversal  lithotripsy 
systems  should  be  investigated  in  the  laboratory  under  finite-amplitude  conditions 
to  evaluate  their  performance.  Time  reversal  mine-hunting  systems  appear  to  be 
unlikely  to  yield  a  real  device  capable  of  destroying  today’s  mines.  This  is  not  due 
to  a  deficiency  in  the  targeting  or  location  of  the  mines,  but  due  to  the  enormous 
acoustic  pressures  that  would  be  required  to  cause  mine  neutralization.  Such  high 
pressures  are  currently  unobtainable  using  time  reversal  arrays. 

For  hyperthermia,  an  original  numerical  simulation  tool  was  designed  and  used  to 
study  the  effects  of  nonlinearity,  inhomogeneity,  and  time-varying  background  cou¬ 
pling  between  the  acoustic  and  the  thermal  fields.  The  temperature  of  the  hot  spot 
was  found  to  be  enhanced  by  nonlinear  generation  of  higher-frequency  harmonics  as 
expected  from  results  found  in  modern  literature  on  the  subject.  Effects  of  tissue 
inhomogeneity  were  studied,  where  an  overall  conservation  of  thermal  deposition  of 
energy  was  observed,  but  an  increase  in  the  peak  local  temperature  was  observed  for 
the  inhomogeneous  media.  Finally,  a  coupling  of  the  acoustic  and  thermal  response 
of  tissue  was  taken  into  consideration.  A  feedback  relation  between  acoustic  and 
temperature  fields  was  obtained,  with  a  resulting  increase  in  the  absorption  coeffi¬ 
cient  and  temperature.  A  lack  of  measured  data  for  tissue  properties  as  a  function 
of  temperature  was  encountered.  This  study  used  the  available  data  for  sound  speed 
and  absorption  coefficient,  but  measurements  of  the  temperature-dependence  of  the 
other  acoustic  and  thermal  properties  of  tissue  would  be  very  helpful  to  future  mod¬ 
elers,  and  could  perhaps  be  combined  with  simulation  tools  to  provide  realtime  FUS 
treatment  planning  on  a  case-specific  basis. 
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Appendix  A 


THE  FDTD  METHOD  AND  OTHER  NUMERICAL 

CONSIDERATIONS 

A.l  Solving  the  Wave  Equation  Using  FDTD 

A  numerical  solution  of  the  wave  equation  may  be  computed  using  the  finite-difference 
time-domain  (FDTD)  method,  which  was  introduced  into  use  in  its  contemporary 
form  by  Yee  [76].  The  FDTD  method  provides  accurate  and  detailed  solutions  down 
to  the  smallest  scales  of  the  problem.  This  strength  of  the  FDTD  method  is  also  its 
greatest  weakness,  as  the  detailed  calculations,  called  full  wave  simulations,  require 
a  large  amount  of  computer  storage  space.  Since  the  wavelengths  in  electromagnetic 
applications  are  large  enough  so  that  most  systems  can  be  represented  by  a  few 
wavelengths,  the  FDTD  method  is  widely  used  in  the  field  of  electromagnetics.  With 
the  development  of  faster  and  larger  computers,  however,  more  practical  problems  in 
acoustics  will  be  solved  using  the  FDTD  method. 

A.  1.1  Approximating  Partial  Derivatives  Using  Finite  Differences 

The  FDTD  method  approximates  the  spatial  and  temporal  partial  derivatives  with 
discrete  differences  obtained  from  Taylor  series  expansions 
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about  some  node  on  the  computational  grid  indexed  in  time  by  n  and  in  space  by  i. 
In  two  spatial  dimensions 

0  =  i(pr«-2p?+pr.)+o(©, 

%  =  -  ‘‘p?-*  +  pr") + 

|^  =  4(p?'^‘-2p?+Pr')  +  C’('5?),  (centered)  (A.2) 

=  4 (2/^”  -  +  4^”"^  -  pT^)  +  ^(<^4 )’  (right-sided) 

W  "  2f  ’’’  ^ 

Note  that  two  versions  of  the  second  derivative  were  given  in  (A.2).  The  reason  for 
this  is  that  in  the  explicit  method  the  future  value  of  pressure  is  isolated  to  one 
side  of  the  equation,  and  is  solved  for  in  terms  of  present  and  past  values  only.  The 
time  derivative  appearing  in  the  D’Alembertian  in  the  wave  equation  was  chosen  to 
contain  the  future  unknown  pressure  term  in  the  FDTD  solution.  The  other  time 
derivatives,  including  the  second  order  derivative  in  the  nonlinear  term  were  all  done 
using  right-sided  finite  differences,  so  they  only  contained  known  quantities.  The 
reason  for  chosing  the  centered  difference  for  the  D’Alembertian  time  derivative  is 
that  the  numerical  solution  behaved  better  with  this  combination. 

Note  that  as  the  order  of  the  derivative  increases,  more  instances  of  the  pressure 
field  are  needed  to  be  held  in  storage.  So  for  example,  the  third  order  time  derivative 
requires  the  current  (n)  pressure  field,  as  well  as  the  pressure  fields  at  previous  time 
steps  (n  —  1), . . .  ,{n  —  5).  This  requirement  is  further  increased  if  the  order  of  the 
differencing  were  higher  than  2.  The  particular  formulae  used  to  approximate  the 
derivatives  can  be  obtained  by  following  one  of  several  formal  methods  [8].  The  exact 
form  will  differ  if  the  differencing  is  to  be  centered  or  right  or  left  handed,  meaning  that 
the  expression  will  contain  indices  around  n  or  to  the  left  or  right  of  n.  Typically,  time 
derivatives  are  formulated  so  that  only  the  known  pressures  from  previous  time  steps 
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axe  required  for  the  calculation,  as  the  explicit  method  steps  forward  and  calculates 
at  each  time  step.  This  technique  is  called  an  explicit  method.  The  first  derivative 
can  be  used  to  calculate  the  difference  expressions  for  higher  order  derivatives  to  the 
same  accuracy  by  applying  the  first  derivative  expression  recursively.  The  spatial 
differencing  is  similar,  but  is  done  using  centered  differences,  and  to  fourth  order 
accuracy,  0{5xY: 

%  ^(-pr+2,i+8pf+i,,.-8pr-ij+p?-2j),  (A.3) 

d^v  1 

si  “  +  +  (A.4) 

The  above  expressions  can  be  found  in  Ferziger  and  Peric  [26]. 

Using  these  centered-difference  expressions  we  can  obtain  the  third  derivative  by 
using  the  first  derivative  expression  as  an  operator  on  the  second  derivative,  yielding 
the  centered  expression  to  fourth  order, 

W  =  (A.6) 

-  248.^jti  +  248*-i  -  158*_2  +  24.^._3  - 

Turning  (A. 5)  into  a  right-sided  finite  difference  expression  can  be  done  if  the  original 
fourth-order  derivatives  are  right-sided,  but  was  not  used  because  it  requires  excessive 
run  time.  One  point  to  make  is  that  the  order  of  the  method  does  not  imply  or 
guarantee  accuracy,  it  only  indicates  how  rapidly  the  truncation  errors  decrease  with 
decreasing  5^  or  5t.  For  this  reason  it  is  possible  for  a  second-order  solution  to  be 
more  accurate  than  a  fourth-order  solution.  While  second-order  time  differencing  was 
used  exclusively  in  this  study,  both  second  and  fourth-order  spatial  differencing  was 
used.  The  highly-steepened  waveforms  in  the  nonlinear  time  reversal  study  behaved 
better  with  fourth  order  differencing  in  space.  Less  numerical  dispersion  (numerical 
propagation  of  different  wavenumbers  at  different  velocities  on  the  grid)  was  observed 
when  the  spatial  grid  spacing  was  made  as  small  as  possible,  and  the  fourth-order 
differencing  was  used. 
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The  need  to  reduce  spatial  grid  size  leads  to  a  requirement  to  reduce  the  time 
grid  spacing  accordingly  for  stability  reasons.  In  Cartesian  coordinates,  the  explicit 
finite-dilference  solution  of  the  lossless  linear  wave  equation  is  stable  only  if 

Cr..Jt  <  .  (A.6) 

This  is  a  statement  that  no  wave  propagation  information  can  travel  across  more  than 
8x  in  a  single  time  step  8t.  The  stability  criterion  becomes  much  more  complicated 
for  the  absorbing  case,  and  simple  Von  Neumann  stability  analysis  is  not  possible  for 
the  nonlinear  wave  equation,  and  is  beyond  the  scope  of  this  study.  However,  it  was 
empirically  observed  that  a  trade  off  between  a  and  8t  exists  in  the  stability  space  for 
the  absorbing  wave  equation.  In  other  words,  increasing  a  necessitates  increasing  8t 
for  stability.  Perhaps  the  only  positive  side-effect  of  the  instability  issue  which  plagues 
explicit  methods  is  that  once  a  code  is  stable,  its  output  is  generally  convergent.  In 
our  studies  we  compared  second  and  fourth-order  spatial  differencing  and  found  very 
good  agreement.  A  full  investigation  of  the  dispersion  and  the  stability  characteristics 
of  the  FDTD  codes  used  is  beyond  the  scope  of  this  study.  Numerical  dispersion  will 
occur  in  all  finite  difference  codes  on  uniform  grids  [64],  but  its  effect  is  minimized 
by  proper  selection  of  the  spatial  grid  size  relative  to  the  dominant  wavelength.  This 
“magic  size”  8t  as  it  is  referred  to  by  Taflove  turns  out  to  be  when  equation  (A.6)  has 
an  equality  sign  for  the  relationship. 

This  study  also  used  implicit  solutions  of  the  wave  equation,  and  the  results  were 
similar  to  the  explicit  results.  However,  the  implicit  solutions  require  more  care  to 
ensure  convergence,  as  large  8t  will  gradually  reduce  accuracy.  By  contrast  increasing 
8t  in  the  explicit  solutions  will  lead  to  instability  and  failure  of  the  computer  code. 

A.2  Modeling  the  Sources 

Both  discrete  axrays  and  continuous  (discretized)  extended  sources  can  be  modeled 
with  the  techniques  used  in  this  study.  To  model  an  array  having  inter-element 
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separation  larger  than  the  spatial  grid  size  a  gap  in  the  placement  of  point  sources  is 
allowed  from  which  no  sound  will  emanate.  To  model  a  continuous  extended  source, 
all  grid  points  lying  on  the  source  are  driven  in  phase.  An  extended  source  described 
by  many  sources  much  smaller  than  a  wavelength  located  along  the  surface  of  the 
simulated  extended  source  can  be  expected  to  collectively  yield  a  wavefront  similar 
to  that  of  the  extended  source  by  Huygens’  construction  [61]. 


(a) 


(b) 


Figure  A.l:  The  layout  of  the  2-D  simulations  in  (a)  Cartesian  coordinates  and  (b) 
polar  cylindrical  coordinates.  The  computational  domain  is  shaded.  The  boundary 
conditions  used  at  the  edges  of  the  domain  are  either  absorbing  boundary  conditions 
(ABC)  or  reflected  boundary  conditions  (RBC). 


A. 2.1  Boundary  Conditions 

Two  types  of  boundary  conditions  are  described  in  this  section.  The  absorbing  bound¬ 
ary  condition  to  prevent  artificial  wave  reflections  from  the  edges  of  computational 
domains,  and  the  reflected  boundary  condition  to  enforce  symmetry  about  a  boundary 
which  lies  on  an  axis  of  symmetry.  Since  the  spherical  bowl  sources  do  not  guarantee 
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that  the  discretized  point  sources  lie  on  exact  grid  positions  in  rectilinear  coordinates, 
the  nearest  grid  location  was  chosen  in  this  case. 


Absorbing  Boundary  Conditions 

Since  numerical  simulations  are  carried  out  on  finite  computational  domains,  the 
issue  of  how  to  handle  waves  reaching  the  outer  edges  of  the  domains  needs  to  be 
addressed.  Untreated,  acoustic  waves  reaching  the  edge  of  a  computational  domain 
would  be  reflected  by  the  artificial  boundaries  of  the  simulation.  In  order  to  avoid 
this  unnatural  behavior  when  modeling  extended  regions  in  space,  it  is  common  to 
use  absorbing  boundary  conditions  (ABC’s)  along  these  edges. 

Several  techniques  are  known  which  prevent  waves  incident  on  domain  boundaries 
from  reflecting.  We  describe  two  different  techniques  here  which  were  used  successfully 
with  the  FDTD  solution  of  the  acoustic  wave  equation.  The  first  involves  placing  an 
absorbing  boundary  layer  along  the  edges  of  the  computational  domain  in  such  a  way 
as  to  minimize  reflected  wave  amplitudes.  This  technique  is  described  in  Kosloff  and 
Kosloff  [50]  for  a  wave  equation  in  2-D, 


2-y^ 


■ Source.  (A. 7) 


This  ABC  requires  padding  the  computational  domain  with  at  least  10  to  20  highly 
absorbing  grid  cells,  on  which  no  useful  calculations  can  be  performed  due  to  the 
artifice  of  the  padding.  Further,  the  slope  of  the  profile  of  7  is  a  compromise  between 
two  evils:  on  one  hand  if  the  profile  is  too  shallow,  a  wider  padded  region  is  necessary 
to  absorb  the  wave.  On  the  other  hand  if  the  absorbing  region  has  a  very  steep  profile 
the  wave  will  be  partially  reflected  by  the  change  in  impedance,  depending  on  the 
wavelength  of  the  wave  compared  to  the  width  of  the  transition  region.  Which  means 
that  the  absorbing  boundary  layer  ABC  results  in  different  reflection  coelficients  for 
different  incident  frequencies,  and  would  have  to  be  optimized  for  best  effect.  A  profile 
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for  7  which  produces  good  results  is 

7  =  17o/cosh^(<7  n), 


(A.8) 


where  Uo  is  a  constant  describing  the  magnitude  of  the  absorption  near  the  edge  of 
the  domain,  and  a  controls  the  slope  of  the  profile,  and  n  is  the  number  of  grid  cells 
over  which  the  transition  takes  place  from  the  nonabsorbing  interior  to  the  edge. 

The  second  ABC  implemented  successfully  in  this  study  uses  a  radiation  condition 
normal  to  the  boundaries  and  was  first  described  for  the  FDTD  method  by  Mur  [60], 
and  is  commonly  known  as  Mur’s  ABC.  This  is  the  ABC  of  choice  for  the  simulations 
presented  in  this  dissertation.  Mur’s  first-order  ABC  uses  the  first-order  radiation 
condition 

dp'  1  dp' 
dx  Co  dt 

along  the  x  =  0  boundary,  and 


dp'  1  dp' 
dx  Co  dt 


=  0  (A.IO) 

X=Xrnax 

boundary.  The  radiation  conditions  ensure  that  waves  normally 


along  the  x  =  x, 

incident  upon  the  boundaries  are  absorbed.  However,  for  obliquely-incident  waves  the 
reflection  coefficient  increases  as  the  direction  of  incidence  departs  from  normal.  For 
grazing  incidence,  the  first-order  Mur  ABC’s  are  very  inefficient  and  a  second-order 
Mur  ABC  may  be  used  [60]. 

Recently,  a  type  of  ABC  known  as  the  perfectly  matched  layer  method  (PML)  has 
been  introduced,  and  is  considered  to  be  the  most  effective  of  all  known  ABC’s  [12]. 
The  method  is  usually  associated  with  spectral  methods,  and  its  implementation  is 
more  complicated  than  the  other  ABC’s  described  above. 


Reflected  Boundary  Conditions 

Another  type  of  boundary  condition  used  in  this  study  is  the  reflected  boundary 
condition  (RBC).  This  is  the  result  of  truncating  the  computational  domain  at  a 
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region  of  symmetry,  where  two  identical  but  opposing  waves  would  meet.  In  our  case 
we  assumed  azimuthal  symmetry  about  the  central  axis  of  bowl  transducers.  Thus 
only  one  half  of  the  2-D  (a;  r)  plane  need  be  computed,  and  the  computational  domain 
is  r  =  [0,r,„aa;],  X  =  [0,Xmax],  as  shown  in  Figure  A.l  is  used  to  obtain  solutions  on 
r  =  [-r^axyVmaxl  X  =  [0,Xmax]-  This  reduces  the  required  storage  space  by  half,  and 
reduces  the  execution  time.  For  such  cases  the  boundary  condition  along  the  r  =  0 
edge  is 


dr 

dr 


=  0, 


=  0. 


(A.ll) 


As  a  natural  result,  the  term  containing  1/r  in  the  Laplacian  for  polar  coordinates  is 
dropped,  since  it  is  multiplied  by  dp'/ dr. 


Other  types  of  boundary  conditions 

The  RBC  is  closely  related  to  the  periodic  boundary  condition,  which  causes  waves 
reaching  and  exiting  one  face  of  the  computational  domain  to  re-enter  the  computa¬ 
tional  domain  from  the  opposite  end.  For  example,  if  computer  memory  is  limited  we 
can  ^Teuse”  available  space  by  allowing  the  wave  to  enter  from  the  left  of  the  domain 
immediately  upon  exiting  from  the  right  side  of  the  domain.  In  both  the  reflected 
and  the  periodic  boundary  conditions,  the  variables  near  the  edges  of  the  domain 
are  replaced  by  values  taken  from  the  adjacent  and  opposite  sides  of  the  boundary 
respectively. 


A. 3  Simulation  Challenges 

Ideally,  computer  simulations  of  acoustic  fields  would  provide  exact  information  re¬ 
garding  all  field  variables  (pressure,  velocity,  etc.)  over  any  size  computational  do¬ 
main  for  any  length  of  time  desired.  Realistically  of  course,  computer  resources  such 
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as  memory  and  execution  time  limit  the  size  of  any  simulation.  Several  simulation 
parameters  determine  the  extent  of  the  maximum  spatial  domain  and  time  span  at¬ 
tainable. 

The  primary  parameter  that  defines  the  size  of  a  simulation  in  acoustics  is  the 
number  of  spatial  dimensions  and  the  number  of  wavelengths  per  dimension.  Until 
recently,  it  was  not  possible  to  simulate  realistic  wave  propagation  problems  in  three 
dimensions  over  ranges  greater  than  a  few  tens  of  wavelengths.  Even  today,  such 
three-dimensional  simulations  are  beyond  the  reach  of  most  engineering  workstations 
because  the  memory  required  to  hold  the  variables  being  calculated  is  proportional 
to  the  total  number  of  spatial  grid  points  used.  So  for  a  n-dimensional  simulation 
having  Nx  wavelengths  on  a  side  in  the  computational  domain,  and  grid  points 
per  wavelength,  the  total  number,  NtotaU  of  data  points  needed  to  hold  one  scalar 
field  quantity  at  any  instant  in  time  is 

Ntotai  =  (A.12) 

where  m  is  an  integer  usually  between  3  and  7,  depending  on  the  highest  order  of 
the  time  derivatives  taken  on  the  unknown  in  the  wave  equation,  and  the  degree 
of  accuracy  used  in  the  time  differencing.  For  a  3-dimensional  simulation,  having  10 
wavelengths  on  a  side  of  a  computational  cube,  each  wavelength  sampled  by  a  modest 
15  grid  points,  Ntotal  is  greater  than  lO’^  when  the  linear  wave  equation  is  solved  to 
second  order  accuracy  in  time. 

For  the  above  reasons,  unless  the  three-dimensionality  significantly  contributes  to 
the  qualitative  yield  of  the  simulation,  two-  and  sometimes  one-dimensional  simula¬ 
tions  are  the  choice  for  most  research  problems.  Another  item  to  note  is  that  while 
N-y  of  15  to  25  is  reasonable  for  most  sinusoidal  waves,  this  number  must  be  increased 
to  adequately  sample  wavefronts  with  steep  gradients,  such  as  shock  waves. 
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